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SECTION  1 


INTRODUCTION 

Two  robustness  issues  are  considered  in  the  research  project.  One  is  the  command/ 
disturbance  input  uncertainties  and  the  other  is  the  plant  uncertainties.  The  controller  is  to 
be  designed  such  that  the  closed-loop  system  remains  stable  for  all  possible  plant 
perturbations  and  the  error  response  is  minimized  for  all  disturbances  and  commands  in  the 
prescribed  set. 

The  prescribed  set  of  disturbances  and  commands  is  modeled  by  the  designer 
according  to  the  real-world  environment  of  the  system.  If  the  set  of  disturbances  and 
commands  under  consideration  can  be  described  by  a  random  process  then  a  stochastic 
design  approach  like  LQG  [1]  or  Wiener-Hopf  [2]  method  can  be  used  to  minimize  the 
mean  value  of  the  error  energy.  On  the  other  hand,  if  we  do  not  have  any  useful  statistical 
information  about  the  disturbances  and  the  commands,  then  a  new  developed  H“ 
optimization  approach  [3-29]  is  recommended.  The  H°°  optimization  approach  is  a  minmax 
design  method  which  minimizes  the  maximal  error  energy  for  all  possible  disturbances  and 
commands  in  the  prescribed  set. 

The  concept  of  H~  optimization  was  initiated  by  G.  Zames  [3]  in  1981.  Since  then, 
this  new  research  area  has  attracted  many  researchers  [4-29].  The  popularity  of  this 
research  area  is  mainly  because  this  area  appears  to  have  a  major  impact  on  future 
engineering  practice. 

tT°  optimization  technique  is  a  powerful  tool  which  not  only  allows  us  to  consider  a 
more  general  class  of  command  and  disturbance  inputs  than  those  considered  by  LQG  and 
W-H  approaches  but  also  has  the  ability  to  handle  the  plant  uncertainties  [10-29].  H~ 
optimization  technique  uses  the  H°°  norm  of  a  transfer  function  matrix  as  the  measure  of  the 
error  response  of  a  system.  This  reflects  the  maximum  error  energy  which  may  occur  in 
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reality.  YJB  controller  parametrization  [2,1 1,30-32]  is  employed  to  characterize  the  set  of 
the  controllers  which  solve  the  stabilization  problem.  Then  among  the  controllers  in  the  set, 
we  will  find  the  one  that  minimizes  the  error  response  subject  to  control-input  and  stability- 
margin  constraints. 

Two  kinds  of  plant  uncertainties  are  considered.  One  is  the  unstructured  uncertainty 
and  the  other  is  the  real  parameter  perturbation.  The  maximum  singular  value  of  the 
complementary  sensitivity  function  is  used  as  a  measure  of  robust  stability  for  the 
unstructured  perturbation  and  the  real  structured  singular  value  (or  the  multivariable 
stability  margin)  is  used  to  measure  the  robust  stability  for  a  system  under  parameter 
perturbations. 

The  control  problem  is  formulated  in  a  general  way.  The  models  of  the  plant 
uncertainties  and  the  reference  and  disturbance  inputs  are  general  enough  to  reflect  the 
situations  found  in  actual  problems  and  the  controller  structure  is  the  most  general  one 
under  linear  time-invariance  which  includes  the  two-degree-of-freedom  controllers. 

Although  our  design  basically  is  a  frequency-domain  approach,  the  computation  is 
not  necessarily  done  in  frequency  domain.  As  a  matter  of  fact,  many  algorithms  developed 
on  the  state- space  framework  have  better  numerical  properties  than  their  counterparts  based 
on  polynomial  matrix  manipulations.  All  the  computations  in  our  design  are  implemented 
on  the  state-space  framework.  Furthermore,  many  pole-zero  cancellations  can  be  done 
theoretically  without  using  any  numerical  minimal  realization  technique.  This  theoretical 
pole- zero  cancellation  technique  [1 1,12,21-24]  can  save  a  lot  of  computing  time  and  avoid 
incomplete  cancellation  due  to  rounding  errors.  It  is  possible  to  perform  all  the  computation 
in  a  numerically  stable,  reliable,  and  efficient  way. 
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Throughout  the  report,  both  of  the  notations 


and  {  A,  B,  C,  D  }  are 


used  for  the  same  purpose  to  represent  a  state-space  realization  of  a  system  whose  transfer 
function  is  C(sI-A)'1B  +  D.  R(s)pxq  is  the  set  of  proper  rational  matrices  with  real 
coefficients.  (RH°°)mxr  is  the  set  of  mxr  proper  rational  matrices  with  real  coefficients 
which  are  analytic  in  the  closed  right  half  plane.  (RL00)1””  is  the  set  of  mxr  proper  rational 
matrices  with  real  coefficients  which  are  analytic  on  the  imaginary  axis. 


Section  2  shows  how  to  formulate  a  control  problem  into  the  standard  H“ 
optimization  problem.  A  mixed-sensitivity  problem  and  a  tracking/disturbance  reduction 
problem  are  employed  to  illustrate  the  formulation  procedure. 

Section  3  reveals  some  useful  properties  of  the  observer-based  controller 
parametrization.  The  poles  of  the  closed-loop  system  with  the  observer-based  controller 
parametrization  are  the  regulator  poles,  the  observer  poles,  and  the  poles  of  the  added 
parameter  matrix.  The  set  of  uncontrollable  and/or  unobservable  controller  poles  is  a  subset 
of  the  regulator  and  the  observer  poles.  The  poles  of  the  closed-loop  system  with  the 
minimal  order  controller  will  include  all  the  poles  of  the  parameter  matrix  and  some  of  the 
regulator  and  the  observer  poles  which  are  not  the  removable  controller  poles.  These 
properties  could  be  useful  in  controller  reduction. 

Section  4  presents  solutions  to  the  standard  H°°  optimization  problem.  In  Sec.  4.1, 
the  observer-based  controller  parametrization  is  used  to  characterize  the  set  of  stabilizing 
controllers  and  to  reduce  the  transfer  matrix  of  interest  as  an  affine  function  of  a  stable 
parameter  matrix.  The  affine  equation  is  then  reduced  further  by  inner-outer  factorizations 
to  a  simple  linear  equation.  In  Sec.4.2,  a  size  reduction  technique  is  used  to  reduce  the 
order  of  the  rational  matrices  involved  in  the  equation.  A  fast  y-iteration  algorithm  for  the 
solution  of  the  two-block  H"  optimization  problem  is  summarized  in  Sec.4.3.  In  Sec.4.4, 
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the  recently  developed  DGKF  approach  [29J  is  employed  to  solve  the  four-block  H“ 
optimization  problem.  Some  numerical  difficulties  we  encountered  will  be  discussed. 

Section  5  presents  our  results  in  the  computation  of  the  real  structured  singular 
value  (or  the  real  multivariable  stability  margin).  In  Sec.5.1,  the  basic  concepts  of  the 
proposed  algorithm  are  briefly  described.  In  section  5.2,  the  mapping  properties  trom  the 
parameter  space  to  the  coefficient  space  and  the  theories  on  which  the  iterative  algorithm  is 
developed  will  be  demonstrated.  The  iterative  algorithm  for  computing  the  real  structured 
singular  value  is  summarized  in  Sec. 5. 3.  Sec.5.4  is  a  summary  of  the  fast  segment  stability 
checking  algorithm  which  is  used  repeatedly  in  the  proposed  iterative  algorithm  in 
computing  the  real  SSV. 

Section  6  is  the  conclusion  and  further  research  suggestions. 
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SECTION  2 


H”  FORMULATION  OF  CONTROL  PROBLEMS 

2.1  Standard  H°°  Optimization  Problem 

Most  control  problems  can  be  formulated  as  the  following  standard  H“  optimization 
problem.  In  the  standard  H~  optimization  problem  formulation,  the  system  representation  is 
rearranged  as  follows. 

z(s)  Gn(s)  G12(s)  v(s)  v(s) 

=  :=  G(s)  (2.1-1) 

_y(s)J  |_g2i(s)  G22(s)J  [u(s)J  Lu(s). 

where  Gn(s)  e  R(s)pxq,  G12(s)  e  R(s)pxm,  G21(s)  e  R^)™*1,  and  G22(s)  e  R(s)rxm 
R(s)pxq  is  the  set  of  proper  rational  matrices  with  real  coefficients.  In  (2.1-1),  z,  y,  v,  and 
u  are  the  controlled  output,  the  measured  output,  the  exogenous  input,  and  the  control  input 
respectively.  The  controlled  input  vector  z  usually  includes  the  error  signal  and  a  weighted 
control  input.  The  exogenous  input  v  contains  the  disturbances,  the  noises,  and  the 
commands.  The  measured  output  vector  y  consists  of  all  the  signals  which  can  be  measured 
and  available  for  feedback.  Through  the  control  input  u,  the  behavior  of  the  system  can  be 
modified.  The  vector  y  will  be  used  as  the  input  to  a  controller  K(s)  and  the  output  of  K(s) 
will  be  connected  to  the  control  input  u.  That  is, 

u(s)  =  K(s)  y(s)  (2.1-2) 

The  standard  H“  optimization  problem  is  the  problem  of  finding  a  proper  controller  K(s) 
such  that  the  closed-loop  system  is  internally  stable  and  IIOII^  is  minimized  where 

IIOIL  =  sup  o  [O(jco)]  (2.1-3) 

0} 

and 

<D(s)  =  G„(s)  +  G12(s)  Q(s)  [  I  -  G22(s)  Q(s)  ]  4  G21(s)  (2.1-3) 
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That  is,  O(s)  is  the  transfer  function  of  the  closed-loop  system  from  v  to  z.  cr[<I>(jct>)]  is  the 
maximum  singular  value  of  O(jco). 

2.2  Formulation  of  Control  Problems  as  the  Standard  H°°  Optimization 
Problem 

As  mentioned  in  the  previous  subsection,  most  of  the  control  problems  can  be 
formulated  as  the  standard  H°°  optimization  problem.  For  the  purpose  of  demonstration, 
two  examples  are  given  in  the  following.  The  first  is  the  mixed-sensitivity  optimization 
problem  and  the  second  is  the  tracking  and  disturbance  reduction  problem. 

2.2.1  A  Mixed-Sensitivity  Optimization  Problem 

Consider  the  system 

y(s)  =  P(s)u(s)  +  v(s)  (2.2- la) 

u(s)  =  K(s)y(s)  (2.2- lb) 

The  nominal  plant  transfer  function  P(s)  is  given  and  the  set  of  disturbances  is  assumed  to 
be  any  signal  vector  with  energy  less  than  or  equal  to  one.  One  of  the  objectives  of 
feedback  control  is  to  find  a  controller  K(s)  such  that  the  closed-loop  system  is  stable  and 
the  maximal  energy  of  the  disturbance  response  is  minimized  subject  to  stability- 
robustness  and  control  input  constraints. 

The  transfer  function  from  v(s)  to  y(s)  is  given  by 

[I-P(s)K(s)]-1  (2.2-2) 

which  is  the  sensitivity  function  of  the  closed-loop  system.  From  [33],  we  have 

il  [I  -  PK]  4  IL  :=  sup  o  [  ( I  -  PK  )4  (jco)  ]  =  sup  {  IlylU  I  HvlU  £  1  }  (2.2-3) 

CD 

That  is,  the  maximal  energy  of  the  disturbance  response  is  equivalent  to  the  H°°  norm  of  [  I 


According  to  Doyle  et.  al.  [34,35],  the  stability  robustness  is  inversely  proportional 
to  the  maximum  singular  value  of  the  complementary  sensitivity  function,  <J[PK(I-PK)' 
1  (jco)].  In  other  words,  the  stability  robustness  is  better  whenever  the  H“  norm  of  the 
weighted  complementary  sensitivity  function  is  smaller.  Usually  the  H“  norms  of  Wj(I- 
PK)1  and  W2PK(I-PK)"!  cannot  be  made  small  at  the  same  time,  i.e.,  if  we  make  one  of 
them  smaller  then  the  other  will  become  larger.  To  have  a  trade-off  between  these  two 
quantities,  Kwakemaak  [14]  formulated  the  mixed-sensitivity  problem  as  the  problem  of 
finding  a  controller  K(s)  which  stabilizes  the  closed-loop  system  and  minimizes  NOIL 
where  <X>  is  given  by 

Wj(I-PK) 1 

<D  = 

-l  (2.2-4) 

W2PK(I-PK) 

W  i  and  W2  are  weighting  matrices  chosen  by  the  designer  according  to  the  real-world 

environment.  If  the  disturbances  most  likely  to  occur  are  low  frequency  disturbances  then 
W^s)  is  chosen  to  be  a  low-pass  filter.  If  the  plant  uncertainty  is  a  function  of  frequency, 
then  W2(s)  should  be  chosen  to  be  frequency-dependent  in  the  same  way. 

In  (2.2-4)  only  the  quantities  related  to  the  disturbance  reduction  and  the  stability 
robustness  are  involved.  The  control-input  constraint  has  not  been  brought  into  the  picture. 
Usually,  the  control-input  energy  is  reduced  when  the  stability  robustness  is  improved  by 
scaling  W^s)  and  W2(s).  If  that  is  not  the  case,  we  can  add  one  term  to  (2.2-4).  That  is, 


O  = 


Wj  (I-PK)  1 
W2PK(I-PK)*1 
W3Ka-PK)4 


(2.2-5) 


where  K(I-PK)'1  is  the  transfer  function  from  v  to  u  and  W3  is  a  weighting  matrix. 
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The  problem  of  finding  a  K(s)  which  stabilizes  the  closed-loop  system  and 
minimizes  IIOII^  can  be  rearranged  into  the  standard  H“  optimization  problem.  Consider  the 

following  system: 


(2.2-6a) 


u  =  Ky  (2.2-6b) 

It  is  easy  to  show  that  the  matrix  G>  defined  by  (2-8)  is  just  the  transfer  function  from  v  to 
[ZjT  z/  z3T]T  of  the  closed-loop  system  (2.2-6).  Comparing  (2.2-6a)  with  (2.1-1),  we 

can  see  that 


(2.2-7) 


2.2.2  A  Tracking/Disturbance  Reduction  Problem 


Fig.  2.2-1  Formulating  a  tracking  problem  as  a  standard  H°°  optimization  problem. 

Consider  the  system  shown  in  Fig.2.2-1.  The  objective  is  to  find  a  controller  K(s) 
such  that  the  closed-loop  system  is  stable  and  the  plant  output  yx  follows  the  command  r  as 
closely  as  possible  under  the  influence  of  the  disturbance  d  and  the  measurement  noise  n. 
Let  the  controlled  output  vector  z  be 
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where  Wj  is  a  weighting  matrix  which  is  chosen  to  constrain  the  control  input.  The 
disturbance,  command,  and  noise  are  assumed  unknown,  but  can  be  described  as 


“d" 

1 - 

£ 

1 _ 

r 

= 

W3 

_n_ 

1 

_ 1 

(2.2-9) 


where  v  is  any  signal  with  bounded  energy.  The  disturbance  d  and  the  command  r  usually 
are  low  frequency  signals  and  hence  W2  and  W3  are  low-pass  filters.  The  weighting  matrix 
W4  is  a  high-pass  filter  since  the  measurement  noise  is  a  high  frequency  signal.  The 
measured  output  y  is 


y  = 


Yj  +  n 
r 


Now,  the  system  can  be  rearranged  as 


’  W2'W3 

P 

Z 

0 

Wl 

y_ 

w2+w4 

P 

1 

JS 

0 

u  =  Ky 


(2.2-10) 


(2.2- 11  a) 


(2.2-1  lb) 


and  the  tracking  problem  is  formulated  as  the  standard  H°°  optimization  problem. 


SECTION  3 


STABILIZATION  AND  OBSERVER-BASED  CONTROLLER 

PARAMETRIZATION 

In  Sec.3.1  the  basic  concept  and  the  historical  developments  of  the  controller 
parametrization  will  be  briefly  described.  Sec.  3.2  will  list  the  previous  results  about  the 
controller  parametrizations  done  by  Youla  et.  al.,  Desoer  et.  al.,  Nett  et.  al.,  and  Doyle  et. 
al..  In  Sec.  3.3,  some  useful  properties  of  the  observer- based  controller  parametrization 
will  be  revealed. 

Throughout  the  section,  the  following  notation  will  be  used.  The  sum,  Ci+B,  of 
two  sets  Cl  with  p  elements  and  B  with  q  elements  is  a  set  which  consists  of  all  elements  of 
Cl  and  B.  Cl+B  has  p+q  elements.  Assume  that  B  C  0,  then  the  difference,  Cl-B,  will 
consist  of  all  the  elements  of  Cl  except  those  in  B.  Cl-B  has  p-q  elements. 

3.1  Basic  Concept  of  Controller  Parametrization 

One  of  the  most  fundamental  requirements  in  control  systems  design  is  to  make  the 
closed-loop  system  internally  stable.  In  addition  to  closed-loop  stability,  usually  the  closed- 
loop  system  is  required  to  meet  some  other  desired  performance  criteria.  Stabilizing 
controller  parametrization  is  important  because  of  the  following  reasons:  (1)  It  provides  the 
full  set  of  the  controllers  which  stabilize  the  closed-loop  system.  (2)  The  full  set  of 
stabilizing  controllers  is  characterized  in  terms  of  a  stable  parameter  matrix  and  the  closed- 
loop  system  is  internally  stable  if  and  only  if  the  parameter  matrix  is  stable.  (3)  The  closed- 
loop  transfer  function  matrix  related  to  the  performance  can  be  written  as  a  simple  affine 
function  of  the  parameter  matrix  and  then  the  control  system  design  problem  becomes  that 
of  finding  a  stable  parameter  matrix  such  that  the  closed-loop  transfer  function  matrix  meets 
the  desired  performance  criteria. 
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The  characterization  of  the  set  of  all  stabilizing  controllers  in  terms  of  a  stable 
parameter  matrix  was  introduced  by  Youla  et.  al.  [2]  in  1976.  Youla's  controller 
characterization  was  developed  based  on  the  fractional  factorizations  over  the  ring  of 
polynomial  matrices.  The  only  drawback  of  Youla's  characterization  is  that  the  stabilizing 
controller  may  not  be  proper.  This  drawback  was  removed  later  by  Desoer  et.  al.  [30]  in 
1980. 

Desoer  et.  al.  [30]  generalized  Youla  et.  al.'s  result  based  on  the  fractional 
factorizations  over  a  general  ring.  The  ring  can  be  chosen  as  the  set  of  proper  stable  rational 
matrices  if  the  given  plant  is  a  linear  time-invariant  system  which  is  represented  by  a 
rational  matrix.  Based  on  the  fractional  factorizations  over  the  ring  of  proper  stable  rational 
matrices  the  set  of  all  proper  stabilizing  controllers  can  be  characterized  in  terms  of  a  proper 
stable  parameter  matrix. 

To  use  Desoer  et.  al.'s  version  of  proper  stabilizing  controller  parametrization,  it  is 
essential  to  compute  the  fractional  factorizations  over  the  ring  of  proper  stable  rational 
matrices.  Nett  et.  al.  [31]  proposed  a  very  convenient  state-space  method  for  this 
computation  in  1984.  The  computation  method  was  developed  based  on  the  observer  and 
regulator  theories. 

Later  in  1984,  Doyle  et.  al.  [11]  showed  that  the  proper  stabilizing  controller 
parametrization  can  be  realized  as  an  observer-based  controller  with  an  added  stable 
parameter  matrix.  The  structure  is  simple  and  easy  to  implement. 

In  Sec.3.3  we  will  show  that  the  observer-based  controller  parametrization  has  the 
following  two  properties:  (1)  The  poles  of  the  closed- loop  system  with  the  observer-based 
controller  parametrization  are  the  regulator  poles,  the  observer  poles,  together  with  the 
poles  of  the  added  stable  parameter  matrix.  (2)  If  the  controller  is  realized  by  a  minimal 
realization,  the  closed-loop  poles  will  include  all  the  poles  of  the  added  stable  parameter 
matrix  and  a  subset  of  the  regulator  and  the  observer  poles. 


The  concept  of  stabilizing  controller  parametrization  is  briefly  described  as  follows. 
Consider  the  block  diagram  in  Fig.  3.1-1  where  v  is  the  exogenous  input  vector  which  may 
consist  of  the  disturbances,  noises,  and  the  commands,  u  is  the  control  input  vector 
through  which  the  behavior  of  the  system  can  be  modified,  z  is  the  controlled  output  vector 
which  is  composed  of  all  the  variables  to  be  controlled,  and  y  is  the  measured  output  vector 
which  consists  of  all  the  measurable  quantities  available  for  feedback.  The  plant  G(s)  is 
given  by 

Gn(s)  G12(s) 

G2i(s)  G22(s) 

The  objective  of  a  typical  control  problem  is  to  find  a  proper  controller  K(s)  which 
stabilizes  the  closed-loop  system  and  the  H"  (or  H2)  norm  of  the  closed-loop  transfer 
function  matrix  <X>(s)  from  v  to  z  is  minimized.  The  first  step  to  solve  the  problem  is  to  find 
the  set  of  all  proper  controllers  which  make  the  closed-loop  system  internally  stable.  Then 
in  the  set  of  all  proper  stabilizing  controllers,  one  will  be  chosen  such  that  the  H~  (or  H2) 
norm  of  <D(s)  is  minimized. 


Fig.  3.1-1  Block  diagram  of  a  typical  control  problem. 


The  stabilizing  controller  parametrization  we  are  interested  in  has  the  following  two 
properties:  (1)  All  the  proper  stabilizing  controllers  can  be  characterized  in  terms  of  a  proper 
stable  parameter  matrix  Q(s)  and  the  closed-loop  system  is  internally  stable  if  and  only  if 
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Q(s)  is  stable.  (2)  The  transfer  function  matrix  <D(s)  is  a  simple  affine  function  of  the 
parameter  matrix  Q(s). 

After  the  stabilizing  controller  parametrization,  the  above  control  problem  becomes 
that  of  finding  a  proper  stable  matrix  Q(s)  such  that  IKDII^  (or  ll<WI2 )  is  minimized.  Property 

(2)  of  the  last  paragraph  is  important  since  it  will  make  the  H°°  (or  H2)  optimization 
problem  easy  to  solve. 

3.2  Youla's  Controller  Parametrization 

The  previous  results  related  to  the  controller  parametrization  will  be  briefly 
reviewed  in  this  subsection.  The  following  theorem  was  originally  developed  by  Youla  et. 
al.  [2]  and  later  modified  by  Desoer  et  al.  [30]. 

Theorem  3.2-1:  [2,  30]  (Youla's  Controller  Parametrization) 

Consider  the  system  in  Fig.  3.1-1.  Assume  that  the  realization  in  (3.1-1)  is  minimal  and 
the  subsystem  G22(s)  is  stabilizable  and  detectable.  Let  M^s),  N2(s),  X2(s),  Y2(s), 
Mj(s),  Nj(s),  Xj(s),  and  Yj(s)  be  proper  stable  rational  matrices  such  that 


I  0 

(3.2-1) 

0  I 

M2(s)-‘N2(s)  =  G22(s)  (3.2-2) 

Then  the  set  of  all  proper  stabilizing  controllers  can  be  described  as 

{  K(s)  I  K(s)  =  [  Mj(s)Q(s)  +  Y2(s)  ]  [  Nj(s)Q(s)  -  X2(s)  ]  A 

with  Q(s)  proper  stable  and  I  Nj(«>)Q(<»)  -  X^**)  |  *0  }  (3.2-3) 


Mjfs)  N2(s)  X2(s)  -Nj(s) 
-Y^s)  Xj(s)  Y2(s)  Mj(s) 


and  the  closed-loop  transfer  function  matrix  ®(s)  from  v  to  z  is  an  affine  function  of  the 
parameter  matrix  Q(s), 
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O(s)  =  Gn(s)  -  G12(s)Y2(s)M2(s)G21(s)  -  G12(s)M1(s)  Q(s)  M^G^s)  (3.2-4) 

To  use  Theorem  3.2-1,  we  need  to  construct  the  proper  stable  rational  matrices  in 
(3.2-1)  and  (3.2-2).  Nett  et.  al.  [31]  proposed  a  convenient  state-space  approach  for  this 
construction.  That  is,  the  following  realizations 


and 


M^s)  N2(s) 
-Y^s)  Xj(s) 


'  a+hc2 

H 

b2+HD22  ■ 

= 

C2 

I 

D22 

-F 

0 

I 

(3.2-5a) 


X2(s) 

Y2(s) 


-N^s) 

M^s) 


a+b2f 

H 

B2 

= 

-(C2+D22F) 

I 

-d22 

F 

0 

I 

(3.2-5b) 


are  proper  stable  and  satisfy  (3.2-1)  and  (3.2-2)  where  F  and  H  can  be  arbitrarily  chosen 
such  that  A+B2F  and  A+HC2  are  stable. 

Doyle  et.  al.  [1 1]  showed  that  if  (3.2-5a)  and  (3.2-5b)  are  used  to  realize  the  proper 
stable  rational  matrices  in  (3.2-1)  and  (3.2-2)  and  let 


rj„M 

J,2(s) 

'  a+b2f+hc2+hd22f 

-H  -(B2+HD22)' 

J(s)  = 

><s> 

J22(s[ 

— 

F 

-(C2+D22F) 

0  -I 

1  °22 

then  the  set  of  proper  stabilizing  controllers  described  in  Theorem  3.2- 1  will  have  a 
structure  as  that  shown  in  Fig.  3.2-1. 


U  «4 


y 


- 1 

h* - 

y 

Q(s) 

with  Q(s)  proper  stable  and  I  -  D22Q(<»)  invertible. 
Fig.  3.2-1  Structure  of  stabilizing  controller  parametrization. 


Replace  the  controller  K(s)  in  Fig.  3.1-1  by  the  structure  of  Fig.  3.2-1,  then  the 
closed-loop  system  can  be  redrawn  as  that  shown  in  Fig.  3.2-2. 


z 


i  ig.  3.2-2  The  closed-loop  system  characterized  in  terms  of  a  parameter  matrix  Q(s). 


with  Q(s)  proper  stable  and  I  -  D^QCoo)  invertible 
Fig.  3.2-3  The  observer-based  controller  parametrization. 


In  Fig.3.2-2,  the  open-loop  transfer  function  matrix  from  u2  to  y,  T22(s),  is  zero. 
Therefore,  the  closed-loop  transfer  function  matrix  from  v  to  z,  i.e.,  O(s),  is  a  simple 
affine  function  of  the  parameter  matrix  Q(s).  That  is. 
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<D(s)  =  Tn(s)  +  T12(s)Q(s)T21(s) 


(3.2-7) 


where  the  realizations  of  Tn(s),  T12(s),  T21(s)  are  given  by 


■  A+B2F 

-b2f 

B, 

Tll«  - 

0 

a+hc2 

VHD21 

Cj+Di2F 

-d12f 

D11 

(3.2-8a) 


T2i(s)  ~ 


■  a+b2f 

V 

_C1+D12F 

°12. 

'a+ho2 

Bl+HD2l' 

C2 

D21 

(3.2-8b) 


(3.2-8c) 


Doyle  et.  al.  [11]  also  pointed  out  that  the  structure  of  the  stabilizing  controller 
parametrization  in  Fig.  3.2-1  can  be  realized  as  an  observer-based  controller  with  an  added 
stable  dynamics  Q(s).  The  realization  is  shown  in  Fig.  3.2-3. 


Note  that  in  Fig.  3.2-3  the  block  diagram  inside  the  dotted-line  box  is  the  well- 
known  full-order  observer-based  controller  [1]. 


3.3  Observer-based  Controller  Parametrization  and  Its  Properties 


In  Fig.  3.1-1,  the  internal  stability  of  the  closed-loop  system  depends  only  on 
G22(s)  and  K(s),  i.e.,  the  interconnected  system  shown  in  Fig.  3.3-1. 


u 


G22(s) 

1  * 

1 

1 

K(s) 

g  | 

Fig.  3.3-1  Equivalent  system  to  Fig.  3.1-1  for  internal  stability. 
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In  this  subsection,  the  controller  K(s)  in  Fig.  3.3-1  is  replaced  by  the  block  diagram 
of  Fig.  3.2-3  which  is  the  observer-based  controller  with  an  added  dynamics  Q(s). 

The  following  theorem  shows  that  the  poles  of  the  closed-loop  system  with  the 

observer-based  controller  parametrization  in  Fig.  3.2-3  are  the  regulator  poles  (the 
eigenvalues  of  A+B2F),  the  observer  poles  (the  eigenvalues  of  A+HC2),  together  with  the 

poles  of  the  parameter  matrix  Q(s).  In  the  design  of  the  observer-based  controller,  F  and  H 
are  chosen  such  that  the  eigenvalues  of  A+B2F  and  A+HC2  are  stable.  Therefore,  the 

closed-loop  system  is  internally  stable  if  and  only  if  the  paramete'  matrix  Q(s)  is  proper 
stable.  The  proof  is  quite  straightforward  and  is  done  completely  in  the  state  space  without 
referring  to  the  derivations  used  by  Youla  et.  al.,  Desoer  et.  al.,  and  Doyle  et.  al.. 

Theorem  3.3-1:  (Observer-based  Controller  Parametrization) 

Consider  the  closed-loop  system  in  Fig.  3.3-1.  Assume  that  G^s)  =  {A,  B2,  C2,  D22) 
with  order  n  is  stabilizable  and  detectable  and  the  controller  K(s)  is  replaced  by  the 
observer-based  controller  with  an  added  m-th  order  dynamics  Q(s)  as  shown  in  Fig. 
3.2-3.  Then  the  set  of  the  closed-loop  poles  is  composed  of  the  n  eigenvalues  of 
A+B2F,  the  n  eigenvalues  of  A+HC2,  and  the  m  poles  of  the  added  dynamics  Q(s).  That 

is,  the  set  of  the  closed-loop  poles  is 


where 


and 


=  J>  +  (p  +  Tp 

closed-loop  regulator  observer  Q(s) 


regulator 


observer 


=  {  n  eigenvalues  of  A+B2F  } 
=  {  n  eigenvalues  of  A+HC, } 


*  {  m  Poles  of  Q(s) ) 


(3.3-la) 

(3.3-lb) 
(3.3- lc) 

(3.3-ld) 


Proof:  See  [65,66] 

It  is  well  known  that  in  the  observer-based  controller  design,  the  closed-loop  poles 
are  the  regulator  poles  (the  eigenvalues  of  A+B2F)  and  the  observer  poles  (the  eigenvalues 
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of  A+HC2)  [1].  Theorem  3.3-1  shows  that  the  above  property  still  remains  when  we  add  a 

dynamics  Q(s)  to  the  observer- based  controller  as  shown  in  Fig.  3.2-3.  The  eigenvalues  of 
A+B2F  and  A+HC2  are  still  parts  of  the  closed-loop  poles  after  we  add  Q(s)  to  the 

controller.  Adding  Q(s)  only  introduces  additional  poles  to  the  closed-loop  system  and  the 

added  closed-loop  poles  are  the  poles  of  Q(s).  If  F  and  H  have  been  chosen  such  that 
A+B2F  and  A+HC2  are  stable,  then  the  closed-loop  system  with  the  observer-based 

controller  parametrization  will  be  internally  stable  if  and  only  if  the  parameter  matrix  Q(s)  is 

proper  stable.  From  Fig.  3.2-3,  it  is  easy  to  see  that  the  controller  K(s)  is  proper  if  Q(s)  is 
proper  and  I  -  D^QC00)  is  invertible. 

With  the  observer-based  controller  parametrization,  the  closed-loop  transfer 
function  matrix  from  v  to  z,  i.e.,  <D(s),  is  a  simple  affine  function  of  the  parameter  matrix 
Q(s).  That  is, 

<D(s)  =  Tn(s)  +  T12(s)  Q(s)  T21(s)  (3.3-2) 

where  Tn(s),  T12(s),  and  T21(s)  are  given  by  (3.2-8).  The  added  dynamics  Q(s)  is  a 
proper  stable  rational  matrix  to  be  chosen  such  that  I-D22Q(«>)  is  invertible  and  <D(s)  has 

some  desired  performance.  No  matter  which  Q(s)  is  to  be  selected,  we  always  have  a  clear 
idea  that  the  closed-loop  poles  will  be  the  eigenvalues  of  A+B2F  and  A+HC2  together 

with  the  poles  of  the  added  dynamics  Q(s)  if  the  controller  is  realized  as  that  shown  in  Fig. 
3.2-3. 


Assume  that  the  orders  of  the  plant  and  the  parameter  matrix  Q(s)  are  n  and  m 
respectively.  If  the  controller  is  realized  as  that  shown  in  Fig.  3.2-3,  then  the  order  of  the 
controller  is  n+m  and  the  closed-loop  system  has  2n+m  poles  described  by  the  set  Pcloscd 
loop  *n  (3.3-1).  The  realization  of  the  controller  in  Fig.  3.2-3  may  not  be  minimal.  Suppose 
it  is  not  and  there  are  r  poles  in  the  controller  either  uncontrollable  or  unobservable,  then  the 
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controller  can  be  realized  by  a  minimal  realization  with  order  n+m-r  and  the  number  of 
closed-loop  poles  will  be  reduced  to  2n+m-r. 

In  the  following  theorem,  we  will  show  that  the  uncontrollable  or  unobservable 
poles  of  the  controller  of  Fig.  3.2-3  must  be  the  eigenvalues  of  A+B2F  or  A+HC2  and  the 

closed-loop  poles  will  always  include  all  m  poles  of  the  added  stable  dynamics  Q(s).  If  r 

pole-zero  cancellations  occur  in  the  controller,  then  the  closed-loop  poles  will  include  m 
poles  of  Q(s),  and  2n-r  eigenvalues  out  of  the  set  Ssregulator+  ^observer  which  was  defined 

in  (3.3-1). 

Theorem  3.3-2: 

Consider  the  closed-loop  system  in  Fig.  3.3-1.  Assume  that  G22(s)  =  {A,  B2,  C2,  D22) 

with  order  n  is  stabilizable  and  detectable  and  the  controller  K(s)  is  replaced  by  the 
observer-based  controller  with  an  added  m-th  order  dynamics  Q(s)  as  shown  in  Fig. 
3.2-3.  A  minimal  realization  of  Q(s)  is  given  as  {A,  B,  C,  D).  Define  ^regulator, 

^observer ^  ^Q(s)  ^  (3.3-lb),  (3.3-lc),  and  (3.3-ld)  respectively  and  let 
^removal  =tthc  controller  poles  which  are  either  uncontrollable  or  unobservable)  (3.3-3) 


F  ,  c 

removal 


F  +  F 

regulator  observer 


(3.3-4) 


and  the  closed-loop  system  with  the  minimal  order  controller  will  have  a  set  of  poles 
described  by  the  following 


F  =  F  +  IF  +F  -  F  1 

closed-loop  with  min.  controller  Q(s)  '  regulator  observer  removal  ’ 


Proof:  See  Appendix  A. 


(3.3-5) 
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SECTION  4 


SOLUTIONS  FOR  THE  STANDARD  H“  OPTIMIZATION  PROBLEMS 


Consider  the  system 


Gn(s) 

G21(s) 


G12(s) 

G22(s) 


v(s) 

v(s) 

:=  G(s) 

_u(s)_ 

_u(s)_ 

(4-la) 


u(s)  =  K(s)  y(s)  (4- lb) 

where  Gn(s)  e  R(s)pxq,  G12(s)  e  R(s)pxm,  G2i(s)  e  IR(s)rxq,  and  G22(s)  6  E(s)rxm. 

Recall  that  the  standard  H°°  optimization  problem  is  the  problem  of  finding  a  proper 
controller  K(s)  such  that  the  closed-loop  system  is  internally  stable  and  IIOII^  is  minimized 

where  d>(s)  is  the  transfer  function  of  the  closed-loop  system  from  v  to  z. 


Let 


'  A 

B1  V 

G(s)  = 

Ci 

D11  D12 

C2 

D21  D22 

(4-2) 


be  a  minimal  realization  of  G(s).  Here 


{  A,  B,  C,  D  ) 


(4-3) 


implies  a  state-space  realization  and  G(s)  =  D  +  C  (si  -  A)'1  B. 

According  to  the  dimensions  of  Gn(s),  G12(s),  G21(s),  and  G22(s),  the  standard 
H~  optimization  problem  has  the  following  cases  to  consider:  (a)  p  >  m,  r  <  q;  (b)  p  <,  m,  r 
<  q;  (c)  p  >  m,  r  £  q;  and  (d)  p  £  m,  r  £  q.  Case  (a)  is  referred  to  as  the  four-block  H“ 
optimization  problem.  Cases  (b)  and  (c)  are  referred  to  as  the  two-block  H~  optimization 
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problem.  Case  (d)  is  referred  to  as  the  one-block  H~  optimization  problem.  The  example  in 
Sec.  2.2. 1  is  a  two-block  problem  and  the  example  in  Sec.  2.2.2  is  a  four-block  problem. 

The  one-block  H**  optimization  problem  was  solved  by  Francis,  Zames,  and  Helton 
[51,  Chang  and  Pearson  16),  Safonov  and  Verma  [7]  and  Glover  [8].  A  so-called  y-iteration 
algorithm  was  originally  proposed  by  Doyle  [  10]  for  solving  the  four-block  and  two-block 
H“  optimization  problems.  The  major  burden  of  the  computation  is  in  the  evaluation  of  the 
optimal  H°°  norm.  Chang,  Banda,  and  McQuade  [20]  proposed  a  fast  y-iteration  algorithm 
for  computing  the  optimal  H°°  norm  of  the  two-block  problem.  For  most  problems  we 
encountered,  the  optimal  two-block  H°°  norm  can  be  accurately  (to  a  double  precision 
accuracy)  obtained  within  three  or  four  iterations.  Once  the  optimal  norm  is  computed,  the 
two- block  problem  can  be  easily  transformed  into  a  one-block  problem  and  then  Glover's 
method  [8|  can  be  used  to  construct  an  optimal  controller. 

Although  Chang  and  Pearson  [17]  and  Chu  and  Doyle  [18]  greatly  simplified  the 
original  Doyle's  four-block  H°°  norm  computation  algorithm,  the  simplified  y-iteration 
algorithms  are  still  not  fast  enough  for  practical  applications.  Recently,  Doyle,  Glover, 
Khargonekar,  and  Francis  [29]  presented  a  celebrating  Riccati-equation  type  solution  to  the 
four-block  H“  optimization  problem.  DGKF's  approach  characterizes  all  possible 
stabilizing  suboptimal  H~  controllers  which  order  is  not  higher  than  that  of  the  plant.  The 
major  computation  involved  is  the  solution  of  two  Riccati  equations.  There  are  still  some 
numerical  difficulties  in  the  solution  of  these  P^cati  equations  when  the  controller  is  close 
to  the  optimal. 

In  the  subsection  Sec. 4. 1 ,  the  observer-based  controller  parametrization  is  used  to 
characterize  the  set  of  stabilizing  controllers  in  terms  of  a  stable  parameter  matrix.  With  the 
controller  parametrization,  the  closed-loop  transfer  function  from  v  to  z  is  an  affine  function 
of  the  parameter  matrix.  After  inner-outer  factorization  manipulations,  the  problem  becomes 


2  1 


(4-4) 


Find  a  §(s)  e  (RH00)"1”  such  that  II  &  II  m  is  minimized 
where  &  is  given  by 


Ru+e> 

r*> 


**22 


(RH0”)"1”  is  the  set  of  mxr  proper  rational  matrices  with  real  coefficients  which  are  analytic 

in  the  closed  right  half  plane.  The  above  problem  is  referred  to  as  the  four-block  IT0 

optimization  problem.  The  two-block  H°°  optimization  problem  is  just  a  special  case  of 
problem  (4-4)  with  either  [R^  R^]  =  0  or  [R{2  R22]  =  0.  When  both  [R21  R22]  =  0 

and  [R12  R22]  =  0,  the  problem  is  one-block  problem.  In  Sec.  4.2,  a  state- space  size 

reduction  method  is  used  to  simplify  the  realizations  of  Rn,  R12,  R21,  and  R^.  A  fast  y- 
iteration  approach  for  the  two-block  H"  optimization  problem  is  summarized  in  Sec.  4.3. 
In  Sec.4.4,  the  numerical  difficulties  in  using  the  DGKF  approach  [29]  to  solve  the  four- 
block  H°°  optimization  problem  will  be  discussed. 


4.1  Controller  Parametrization  and  Inner-Outer  Factorization 

Suppose  the  realization  {  A,  B2,  Cj,  D22  }  is  stabilizable  and  detectable,  then  by 
using  Doyle's  version  [11,32]  of  YJB  controller  parametrization  [2,30],  the  set  of 
stabilizing  controllers  can  be  characterized  by  the  structure  shown  in  Fig.3.2-3. 

In  Fig.3.2-3,  the  regulator  gain  matrix  F  and  the  observer  gain  matrix  H  are  chosen 
such  that  A+B2F  and  A+HC^  have  no  eigenvalues  in  the  closed  right  half  of  s-plane.  With 
the  controller  K(s)  parameterized  as  above,  the  closed-loop  transfer  function  matrix  from  v 
to  z,  i.e.,  O(s),  is  a  simple  affine  function  of  the  parameter  matrix  Q(s).  That  is. 


<D(s)  =  Tu(s)  +  T12(s)  Q(s)  T2,(s)  (4.1-1) 

where  the  realizations  of  Tn(s),  Tj2(s),  T21(s)  are  proper  stable  rational  matrices  given 
by  (3.2-8). 
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If  Tj2(s)  and  T21(s)  are  thin  and  fat  respectively,  i.e.,  p  >  m  and  r  <  q,  then  they 
can  be  decomposed  as 

T12(s)  =  T.(s)T(s) ,  T2,(s)  =  t  (s)t(s)  (4.1-2) 

where  T.  e  (RH~)pxm  is  inner,  T.  e  (RH“)rxq  is  co-inner,  Tq  e  (RH~)mxm  and  fQ  e 
(RH00)™  are  outer  matrices.  Furthermore,  we  can  find  a  T±(s)  e  (RH°°)px(p'm)  such  that 
[  T.  T±  ]  is  square  and  inner.  Similarly,  T±(s)  e  (RH°°)(q'r)xq  can  be  found  such  that  [  tf 

-T  T  . 

T±  ]  is  square  and  co-inner.  Now,  (4. 1-1)  can  be  rewritten  as 


®  -  T„  +  [Ti  TJ 


(}  0 

0  0 


here  0  =  TqQ  Tq  .  Define  &  =  [  T.  T±  ]*  O  [  T*  T*  ],  then  we  have 


T?T,,T*  +  Q 

i  11  t  ^ 

£ 

T^T  T* 
i  ln1x 

"  R11  +  ^ 

R12 

T*T..f. 

X  11  l 

T*T  T* 

iXlll1X_ 

Rji 

**22. 

Note  that  II  &  II  m  =  II O II  m  and  therefore  the  problem  becomes  to: 

Find  a  Q(s)  e  (RH”)"1”  such  that  II  &  II  ^  is  minimized 
where  &  is  given  by  (4.1-4). 

Problem  (4.1-5)  is  referred  to  as  the  four-block  H°°  optimization  problem. 


(4.1-3) 


(4.1-4) 


(4.1-5) 


The  two-block  H~  optimization  problem  is  just  a  special  case  of  problem  (4.1-5) 
and  the  corresponding  &(s)  is  a  special  case  of  (4. 1-4)  with  either  [R21  R^]  =  0  or  [R{2 

R^2]  =  0.  When  both  [R21  R^]  =  0  and  [R^2  R22]  =  0,  the  problem  is  one-block 


problem. 

In  (4.1-4),  R(s)  is  given  by 
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R(s)  = 


(4.1-6) 


Rh(s)  Rj2(s) 

=  R.  (s)  RB(s) 

Rj,(»)  R22<S>J  L 


The  state-space  realizations  of  RL(s)  and  Rr(s)  are  given  by  [12] 


-1/2 

is  the  orthogonal  complement  of  D12Rq  so  that 
["  DJ2Rp/2  D±j  is  square  and  orthogonal. 


(4.1 -7a) 

(4.1 -7b) 

(4.1 -8a) 
(4.1 -8b) 

(4.1-8c) 

(4.1-8d) 

(4.1-8e) 
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and 


K  =  -(B,dJ,+YcJ)R^ 


RD  ~  D21D21 


Y  =  Ric 


<A-BiDIiRSC2>T 


BiDIDX 


-CIRDC2 


(A-B15IRDC2> 


D±  is  the  orthogonal  complement  of  RD  D21  so  that 


r  1/2  n 

rd  °21 


D, 


is  square  and  orthogonal. 


(R^KRd2)7  =  rd 


(4.1 -9a) 
(4.1 -9b) 


(4.1 -9c) 


(4.1-9d) 


(4.1-9e) 


In  (4.1-8c)  and  (4.1-9c),  Ric[H]  denotes  the  solution  of  the  Riccati  equation  corresponding 
to  the  Hamiltonian  matrix  H. 
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4.2  Size  Reduction  In  Four-Block  H°°  Formulation 


Neither  RL(s)  nor  Rr(s)  in  (4.1-7)  is  a  minimal  realization.  Limebeer  and  Hung 
[22]  derived  a  minimal  realization  of  Rn(s)  in  the  one-block  case  ( the  case  with  R12(s), 
R21(s)  and  R^s)  being  zero  matrices )  for  finding  a  bound  on  the  McMillan  degree  of  the 
H“  optimal  controller.  This  minimal  realization  process  is  not  numerically  reliable  although 
it  is  adequate  for  their  theoretical  purpose.  Recently,  Postlethwaite,  Gu,  and  O’Young  [23] 

unveiled  some  useful  properties  on  the  solution  of  the  algebraic  Riccati  equation  and  used 
them  to  develop  a  numerically  reliable  algorithm  for  the  minimal  realization  of  Rn(s)  in  the 

one-block  case. 

Partitioning  RL(s)  and  Rr(s)  as 


R,(s)  = 


RL1(s> 

RL2<S>  J 


.  R.(s)  =  [RR1(s>  RR2<S>  ] 


then  we  have 


Rn(s)  =  RL1(s)  Rr1(s)  ,  R12(s)  =  RL1(s)  RR2(s) 


R2i(s)  =  RL2(s)  RRi(s)  ,  R22(s)  =  RL2(s)  RR2(s) 

where 

- 

Rl2(s)  = 


-(A+B2F)T 

(Cj+D12F)T  -XH 

-<b2Rd1/2)t 

<d12Rd1/2)t  0 

-(A+B,F)' 


(Cj+D^F)  -XH 


dJcx! 


D 


0 


(4.2- la) 


(4.2- lb) 


(4.2-  lc) 


(4.2- Id) 
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(4.2-  le) 


-(A+HC2)T 

-(Rd^ 

Rri(s)  = 

Cj  Y+Dj  j  (B !  +HD2  j  )T 

d„(Rd1/2d21)t 

0 

pW 

kd 

RR2(s)  = 

The  recently  discovered  properties  on  the  solution  of  the  algebraic  Riccati  equation  by 
Postlethwaite  et.  al.  [23]  and  the  pole-zero  cancellation  technique  by  the  similarity 
transformation  are  used  for  constructing  the  minimal  realizations  of  RL1(s),  R^Cs),  Rri(s), 
and  RR2(s).  Only  orthogonal  transformations  are  involved  in  the  computation,  so  the 

algorithm  is  numerically  reliable.  The  new  result  of  Postlethwaite  et  al.  is  listed  as  follows. 

Theorem  4.2-1:  [23] 

Consider  the  algebraic  Riccad  equation, 

AtX  +  XA  -  XBBtX  +  CTC  =  0  (4.2-2a) 

There  exists  an  orthogonal  matrix  U  =  [  ^2  ]  such  that 

Ai  Aj 

UTAU  =  (4.2-2b) 

L  0 

t  n 
U  B 

1  ,  CU  =  [0  CU2]  (4.2-2C) 

u2b 

and 
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t  Too 
uxu  =  [o  X, 

The  eigenvalues  of  Aj  are  unobservable  but  stable.  That  is,  the  set  of  eigenvalues  of  A3 
includes  all  the  observable  and  all  the  unstable  eigenvalues  of  A.  The  rank  of  X2,  and 
therefore  the  rank  of  X,  equals  the  dimension  of  Ay  X2  >  0  satisfies  the  following 

algebraic  Riccati  equation  (ARE), 

A^X2  +  X2A3  -  X2U2BBTU2X2  +  U2CTCU2  =  0  (4.2-2e) 

The  algorithm  we  developed  is  different  from  that  of  Postlethwaite  et.  al..  The 
advantage  of  the  decomposition  R(s)  =  Rl(s)Rr(s)  is  taken  to  simplify  the  size  reduction 
process  for  Rn(s),  R12(s),  R21(s)  and  R^s).  Although  we  do  not  try  to  construct 
minimal  realizations  for  R^s),  we  do  have  minimal  realizations  for  Ru(s)  and  RRj(s), 

i,j=l,2,  and  there  is  no  possible  mathematically  identical  pole- zero  cancellation  between 
Ru(s)  and  RRj(s). 

In  the  following,  we  will  derive  minimal  state- space  realizations  for  RL1(s)  and 
RL2(s).  Minimal  realizations  of  RR1(s)  and  RR2(s)  can  be  obtained  in  a  similar  way.  First  of 
all,  let  us  consider  RL1(s). 

Minimal  Realization  of  RL1(s) 

The  first  step  is  to  apply  Theorem  4.2-1  to  the  solution  of  the  algebraic  Riccati 
equation  (4.1 -8c).  There  exists  an  orthogonal  matrix 

U  =  [Uj  U2]  (4.2-3a) 

such  that 

UT(A  -  B2RqD[2C.)U  =  1  12  (4.2-3b) 

2  D  12  1  |_  0  A2 
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where  X2  is  the  positive  definite  solution  of  the  following  reduced-order  algebraic  Riccati 
equation 

A^X2  +  X2A2-X2L2Ljx2  +  W2W2  =  0  (4.2-3e) 

Applying  the  similarity  transformation  (UT,  U)  to  the  realization  of  RL1(s)  in  (4.2- lc)  and 
then  using  (4.1-8)  and  (4.2-3),  we  have 


-Al 

0 

0 

0 

Ru(s)  = 

t 

> 

to 

t 

r 
►— » 

S'* 

X 

to 

-(A2-L2l^X2)T 

(DjW2-D12RoB2U2X2)T 

-XjU^H 

< 

A 

(Di2Rd/2)T 

0 

(4.2-4) 

The  size  of  the  state-space  realization  of  Ru(s)  has  been  reduced  from  the  dimension  of  A 
to  that  of  A2  and  the  realization  in  (4.2-4)  is  controllable  according  to  the  following. 

Remark  4.2-1: 

The  state-space  realization  of  RL1(s)  in  (4.2-4)  is  controllable. 

Proof:  See  Appendix  A. 
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Remark  4.2-2: 


If  R^s)  =  0,  i.e.,  D±  =  0,  then  it  is  easy  to  prove  that  the  pair  {L2T,  (A2-L2L2TX2)T1  is 
observable.  In  this  case,  the  realization  of  RL1(s)  in  (4.2-4)  is  already  a  minimal 
realization. 


Now  the  realization  of  Ru(s)  in  (4.2-4)  is  controllable  but  not  necessarily 
observable  if  RL2(S)  *  0-  By  using  the  Van  Dooren's  version  [36]  of  Rosenbrock's 
algorithm  [37],  we  can  find  an  orthogonal  matrix 


V=[V,  V2] 

such  that 

-vt(vl2lJx/v 

-v](A2-L2LT2x2)Tvl  0 

.-^(V^xA  -v^-L^x/v, 

and 

-  L^V  =  [  -  L^V,  o] 

and  the  pair 

{ .  -VX-L^X/V,  } 


(4.2-5a) 


(4.2-5b) 


(4.2-5c) 


(4.2-5d) 


is  observable. 

Applying  the  similarity  transformation  (VT,  V)  to  the  realization  of  RL1(s)  in  (4.2-4) 
and  eliminating  the  unobservable  part,  we  have  the  following  minimal  realization 
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Minimal  Realization  of  RL2(s) 

The  first  step  is  just  similar  to  that  of  the  minimal  realization  of  RL1(s).  Applying  the 
similarity  transformation  (UT,  U)  to  the  realization  of  R^s)  in  (4.2- Id),  and  then  using 
(4.1-8)  and  (4.2-3),  we  have 


For  the  same  reason  of  Remark  4.2-1,  the  realization  of  R^fs)  in  (4.2-7)  is  controllable. 
Next,  we  can  find  an  orthogonal  matrix 

S  =  [Sj  s2]  (4.2-8a) 

such  that 

-S^Aj-L^X^S 

-sT(V^l£x2)tSi  0 

=  (4.2-8b) 

■S X  -  -sJ(A2-L!^x2)Ts2_ 

and 
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W2X2  S  - 

WjXj's,  0' 

(4.2-8c) 

and  the  pair 

(  WjXjS,  , 

-Sfo  -  L2L2X2)TS1  } 

(4.2-8d) 

is  observable. 

Applying  the  similarity  transformation  (ST,  S)  to  the  realization  of  R^fs)  in  (4.2-7) 
and  eliminating  the  unobservable  part,  we  have  the  following  minimal  realization: 


rl2(s) 


'-s](A2-L2LT2x2)Tsl 

sI<Diwa-Di2RDBIuax/  -S>2U> 

w2x21s1 

0 

HaH 

(4.2-9) 


Minimal  Realization  of  Rr1(s) 

The  procedure  is  similar  to  that  of  Ru(s).  Theorem  4.2-1  will  be  applied  to  the 
solution  of  the  algebraic  Riccati  equation  (4.1 -9c).  There  exists  an  orthogonal  matrix 


u  =  [u,  U4] 


such  that 


uV-BjD^c/u  = 


0 


A 

4  *34 


(4.2- 10a) 


(4.2- 10b) 


UT(R^C2)T  = 


L3 


.  6  b[u  =  [°  w4] 


(4.2-lOc) 


and 


utyu  = 


0  0 

o  Y, 


(4.2- lOd) 


where  Y2  is  positive  definite  solution  of  the  following  reduced-order  algebraic  Riccati 
equation 


32 


(4.2- lOe) 


A4  Y2  +  Y2  A4  -  Y2L4L4  Y2  +  W4  W4  =  0 

Applying  the  similarity  transformation  (UT,  U)  to  the  realization  of  Rrj(s)  in  (4.2-le)  and 
then  using  (4.1-9)  and  (4.2-10),  we  have 


_A3 

-Aj4+  l,l>2 

0 

-a4  +  l4lJy2 

-L4 

0 

# 

d>d%>T 

0 

0 

-a4  +  l4lJy2 

» 

# 

Dh(Rd1/2D21)T 

0 

-1/2 

kd 

• 

where 

#  -  c1u4y2  +  du<d*w4-d*r;Jc2u4y1> 


(4.2- 11a) 


(4.2-llb) 


The  size  of  the  state-space  realization  of  Rri(s)  has  been  reduced  from  the  dimension  of  A 
to  that  of  A4  although  the  realization  in  (4.2-11)  in  general  is  neither  controllable  nor 
observable. 


Remark  4.2-3: 

If  RR2(s)  =  0,  i.e.,  D±=  0,  then  it  is  easy  to  prove  that  the  pair 

{(a4-l4lJy2),l4) 

is  controllable.  In  this  case,  the  realization  of  Rri(s)  in  (4.2-1 1)  is  controllable. 
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There  exists  an  orthogonal  matrix 


v=[v3  V4] 


such  that 


-vV-ljJyjv 


k4  4  4  2' 

-vJ(A4-L4L>2)Y3 


0 


-v>4  -  L4L>2)V3  -v>4  -  L4lJy2)V4 


and 


■#■ 

r#v3  oi 

0 

V  = 

o 

o 

_ 1 

and  the  pair 


(#v3.  -vJ(A4  -  L4L^Y2)V3  ) 


is  observable. 


(4.2- 12a) 


(4.2- 12b) 


(4.2- 12c) 


(4.2-1 2d) 


-T  — 

Applying  the  similarity  transformation  (V  ,  V)  to  the  realization  of  Rri(s)  in 
(4.2- 11)  and  eliminating  the  unobservable  part,  we  have  the  following  realization  which  is 
observable 


1 

« 

1 

*r 

T 

-V  L 

V3L4 

Rr1(s)  = 

#V3 

D..»dX>T 

1 - 

o 

C 

(4.2-13) 


If  Rr2(s)  =  0,  then  from  Remark  4.2-3  the  realization  of  Rri(s)  in  (4.2-13)  is  controllable 
and  therefore  is  minimal.  Otherwise,  we  can  find  an  orthogonal  matrix 

S  =  [S3  S4]  (4. 2- 14a) 

such  that 
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STv[(A4-L4lJV2)V3S 

-sX(a4-l4lJy2)v3S3  -sJvJ(a4-l4lJy2)v3s4 
o  -sJv3(a4-l4lJy2)v3s4 


-slvl 

— t  t  3  3  4 

-SV'L-  0 


(4.2-145) 


(4.2-  14c) 


and  the  pair 


{  -S»4-L4L>2)V3S3  ,  -s^l4  } 


(4.2- 14d) 


is  controllable. 


-T  — 

Applying  the  similarity  transformation  (S  ,  S)  to  the  realization  of  RR1(s)  in 
(4.2-13)  and  eliminating  the  uncontrollable  part,  we  have  the  following  minimal  realization: 


-s3  v3  ( a4  -  l4lJy2  ) v3s3  -s3  v3  l4 


RH1(S)  = 


#v3s3 


d„(Rd/2d21)T 


(4.2-15) 
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Minimal  Realization  of  RR2(s) 

Applying  the  similarity  transformations  (UT,  U)  and  then  (VT,  V)  to  the  realization  of 
RR2(s)  in  (4.2- If)  and  eliminating  the  unobservable  part,  we  have  the  following  realization 
which  is  observable: 

RR2(s)  = 

Furthermore,  we  can  find  an  orthogonal  matrix 

i=[T3  T4]  (4.2- 17a) 

such  that 


(4.2- 17b) 


(4.2- 17c) 


and  the  pair 


{  .tX<a4-l4l*y2)v3t3,  ) 


(4.2- 17d) 


is  controllable. 


Applying  the  similarity  transformation  (T,  T)  to  the  realization  of  RR2(s)  in 
(4.2-16)  and  eliminating  the  uncontrollable  part,  we  have  the  following  minimal  realization: 
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4.3  Fast  Y*'terat*on  Algorithm  for  the  Solution  of  the  Two-block  H“ 
Optimization  Problem 

4.3.1  Basic  Concept  of  Fast  y-iteration 


In  (4.1-4),  if 


&(s)  = 


Rn(s)  +  0(s) 
R^s) 


(4.3-1) 


then  we  have  the  two-block  H~  optimization  problem  described  as  foUows, 

Find  a  0(s)  e  (RH00)"1”  such  that  II&  11^  is  minimized 

where  Misgiven  by  (4.3-1).  (4.3-2) 

The  two-block  H~  optimization  problem  arises  in  the  optimal  disturbance  attenuation  with 
control  weighting,  the  minimization  of  a  weighted  sum  involving  both  the  sensitivity 
function  and  its  complement,  and  any  control  problems  with  more  controlled  outputs  than 

control  inputs.  The  two-block  H°°  optimization  problem  usually  is  solved  in  two  steps.  The 
first  step  is  the  computation  of  the  optimal  H°°  norm  ll&ll  ^  which  is  the  most  time- 

consuming  job.  The  second  step  is  the  construction  of  an  optimal  Q(s)  and  then  an  optimal 
controller. 

In  this  section,  a  very  fast  iterative  algorithm  was  developed  for  the  computation  of 
the  optimal  two-block  H°°  norm,  i.e.,  the  computation  of 


inf 

f)e  (RH“)m” 


R.i+Q 

.  . 


(4.3-3) 


Rn(s)  and  R2i(s)  are  in  (RL°°)mxI  and  (RL~)qxr  respectively.  (RL*0)"1”  is  the  set  of  proper 
rational  matrices  with  real  coefficients  which  are  analytic  on  the  imaginary  axis. 
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In  the  two-block  H°°  optimization  problem,  the  most  computationally  demanding 
work  is  the  computation  of  the  optimal  H°°  norm  y .  Define 

=  be  ll[R,.w+<fc>]M(‘-r>  ’IL  <4-3-4*) 

where 

T2!  -  R2i(-s)T  R2i(s)  =  M(-s,y)T  M(s,y)  (4.3-4b) 

Then  the  problem  of  computing  yo  can  be  considered  as  that  of  finding  a  y  such  that  p(y) 
equals  to  1. 

In  Y-iteration  approach,  an  initial  guess  of  the  optimal  norm  y0,  say  y,  is  made  and 
the  computations  of  the  spectral  factorization  (4.3-4b)  and  the  one-block  optimal  IT*  norm 
computation  (4.3-4a)  are  performed  to  determine  what  the  next  guess  shall  be.  The  next 
guess  of  yis  determined  by  the  value  of  |i.(y).  This  guess  and  computation  loop  is  repeated 
until  p(y)  =  1.  Numerically,  this  computation  ends  when  lp(y)  -  II  <  e,  where  e  is  a  very 
small  number.  This  search  scheme  works  all  right  if  lp(y)  -  II  <  e  implies  that  ly  -  yj  <  e, 

i.e.,  the  absolute  value  of  the  slope  of  p(y)  at  y  is  large  enough.  On  the  other  hand,  if  the 
absolute  value  of  the  slope  is  small,  we  may  have  I  y  -  yc  I  »  e.  In  other  words,  the 
accuracy  of  y  may  be  very  poor  although  the  accuracy  of  p(y)  is  good. 

A  simple  bisection  search  scheme  can  be  used.  Starting  from  a  lower  bound  yL  and 
an  upper  bound  y^,  one  can  define  y  =  (yL  +  yu)/2  and  evaluate  ji(y).  If  p(y)  >  1,  then  yL  is 
updated  by  y.  Otherwise,  yu  is  replaced  by  y.  This  process  is  repeated  until  the  gap  between 
yL  and  Yy  is  small  enough.  This  search  scheme  is  straightforward,  but  the  convergence  rate 
is  slow.  Say,  the  initial  gap  between  yL  and  Yy  is  1.  To  reduce  the  gap  to  be  less  than  10'15, 
the  number  of  iterations  required  is 

n  >  15/(log102)  =  49.8 
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In  Sec.4.3.3  a  very  fast  iterative  algorithm  is  developed  for  the  computation  of  the 
optimal  H°°  norm  yo  in  (4.3-3).  To  achieve  this  objective,  three  tasks  have  been  done.  The 

first  task  is  to  find  better  initial  lower  and  upper  bounds  of  y .  Here  "initial  lower  and  upper 

bounds"  means  the  lower  and  upper  bounds  which  can  be  easily  computed  without  doing 
any  spectral  factorization.  The  second  task  is  to  reduce  the  computation  amount  in  each 
iteration.  The  last,  but  the  most  important  task  is  to  find  a  fast  search  scheme  by  which  the 
number  of  iterations  is  greatly  reduced. 

In  Sec.4.3.2  we  present  an  easy  way  of  computing  a  new  initial  lower  bound  and 
showed  some  useful  properties  of  |i.(y)  which  had  just  been  discovered.  Together  with  the 
well-known  properties  of  p(y),  these  discoveries  are  of  great  help  in  developing  a 
sophisticated  and  fast  search  scheme  for  y  .  The  initial  lower  and  upper  bounds  are  denoted 
by  yt  and  y2  respectively. 

Before  entering  the  iteration  loop,  we  compute  p.(Yj),  and  p(y2).  With  the 
knowledge  of  yv  P-fy),  y2,  and  p(y2),  we  can  generate  three  new  upper-bound  candidates 

and  two  new  lower-bound  candidates  very  easily.  A  new  greatest  lower  bound,  denoted  by 
y3,  and  a  new  least  upper  bound,  denoted  by  y4,  can  be  found. 

Now  we  are  starting  the  iteration  loop  with  the  input  data:  yr  pty),  y2,  p(y2),  yy 
and  y4.  Let  y  =  (y3  +  y4)/2,  and  evaluate  p(y).  If  p(y)  is  exactly  equal  to  1,  then  y  is  the 
optimal  norm  yo.  Otherwise,  p.(y)  may  be  greater  or  less  than  1.  In  either  case,  we  can 

easily  generate  two  new  lower-bound  candidates  and  two  new  upper-bound  candidates  by 
using  the  information  of  y,  n(y),  yr  nty),  y2,  and  n(y2).  The  bounds  y3,  and  y4  can  be 

updated  by  the  new  greatest  lower  bound  and  the  new  least  upper  bound  respectively. 
Depending  on  the  sign  of  p(y)-l,  either  the  pair  {yr  p(y,)}  or  the  pair  {y2,  p(y2)}  will  be 

updated  by  {y,  |i(y) } .  The  iteration  loop  is  repeated  until  the  difference  y4-  y3  is  small 

enough. 
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The  proposed  fast  7-iteration  search  scheme  is  much  faster  than  any  other  existing 
y-iteration  search  scheme.  We  only  require  one  evaluation  of  p(y)  per  iteration.  In  each 
iteration,  two  new  lower-bound  and  two  new  upper-bound  candidates  are  generated.  In 
sevcial  exiuiiples  we  have  encountered,  the  gap  of  lower  and  upper  bounds  can  be  reduced 
from  1  to  10'15  in  only  2  or  3  iterations. 

4.3.2  Bounds  and  Properties  of  the  Function  ^i(y) 

For  notational  simplicity,  we  will  omit  the  indeterminates  s  or  jeo  when  there  is  no 
confusion.  A*  is  the  conjugate  transpose  of  A  if  A  is  a  constant  matrix.  The  same  notation 
A*  is  also  used  as  a  short-hand  notation  for  AT(-s)  when  A  is  a  rational  matrix. 

First  of  all,  we  need  to  find  the  initial  lower  and  upper  bounds  of  y  .  The  initial 

bounds  are  those  which  can  be  easily  obtained  without  doing  any  heavy  computation  like 
spectral  factorization  or  Hankel-Toeplitz  approximation.  The  initial  bounds  which  appeared 
in  the  literatures  [17,18,67]  are  listed  as  follows.  The  initial  upper  bounds  are 


Yui 


‘'11 

*21] 


£  Yn 


(4.3-5a) 


*  Y0 


(4.3-5b) 


where 


and 


a  = 


inf 

(RH° 


>ynxr 


II  Rn+  Oil, 


(4.3-5c) 


b  =  II  R2i  II, 


(4.3-5d) 


The  initial  lower  bounds  are 
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Yli  -  b  £  Yo 


(4.3-6a) 


Yl2  :=  a  ^  Yo 


(4.3-6b) 


.  inf 

§e  (RH00)1"” 

fye  (RH00)^ 


ftl 

61J 


^  Vo 


(4.3-6c) 


where  a  and  b  are  given  in  (4.3-5c)  and  (4.3-5d)  respectively.  The  lower  bound  YL3  is 
greater  or  equal  to  y^. 


Now  we  have  the  best  pair  of  initial  lower  and  upper  bounds  as  follows. 

YL  =  max  {  YL1,  Yl3  )  (4.3-7a) 

and 

Yu  =  min{  Ym,YU2}  (4.3-7b) 

To  design  a  fast  and  effective  iterative  algorithm  for  the  computation  of  yo.  it  is 
important  to  have  enough  knowledge  about  the  properties  of  the  function  p(Y)  which  is 
defined  in  (4.3-4).  Chu,  Doyle  and  Lee  [18],  and  Wang  and  Pearson  [68]  have  made  great 
contributions  in  this.  Their  results  are  listed  in  the  following. 


Theorem  4.3-1:  [18] 

The  function  p(y)  is  a  continuous,  convex,  and  strictly  monotonically  decreasing 
function  of  y. 

Theorem  4.3-2:  [18] 

Define 


f(a,Y)  =  H(a)  j 

'a2-b2/vVV 

(4.3-8) 

for  some  a  >  b,  where  b  =  II  R2i  II.,.  Then 

i)  M-(Y)  <  f(a,Y)  if  Y<«. 

(4.3-9a) 

ii)  p(Y)  =  f(a,Y)  if  Y  =  «- 

(4.3-9b) 

iii)  p(Y)  >  f(a,Y)  if  Y>«- 

(4.3-9c) 
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Theorem  4.3-3:  [68] 
If  |i(Y)  >  1,  then 


Yo  <  M-(Y)  Y 


(4.3-10) 


Theorem  4.3-4:  [68] 

If  (i(y)  <  1,  then 

Y0  <  J  (j,2(y)  y2  +  ( 1  -  |i2(y)  ]  b2  (4.3-11) 

where  b  =  II  R21 11,^. 

The  bound  in  Theorem  4.3-4  can  also  be  derived  from  Theorem  4.3-2.  However, 
Wang  and  Pearson  found  this  bound  before  Chu,  Doyle  and  Lee  obtained  Theorem  4.3-2. 
Chu,  Doyle  and  Lee  never  used  this  bound  in  their  algorithm.  From  Theorem  4.3-2,  it  is 
easy  to  have  the  following  corollary. 

Corollary  4.3-5: 

If  |J.(y)  >  1,  then 

Y°  >  J M-2(Y)  Y2  +  [  1  -  H2(Y)  1  b2  (4.3-12) 

where  b  =  II  R2i  IK 

Proof  of  Corollary  4.3-5: 

M-(Y)  >1  =>  y<yo  =>  from  Theorem  4.3-2  (i)  we  have 

H2(Y)  (Y2  -  b2)  <  Y02  -  b2  =>  (4.3-12). 

The  following  theorem  is  a  modified  version  of  Theorem  4.3-3.  The  proof  is 
similar  to  that  of  Theorem  4.3-3  and  therefore  is  omitted. 

Theorem  4.3-6: 

If  M-(Y)  >  1*  then  Y0  ^  yj p,2(y)  ^  +  [  1  -  p,2(y)  ]  c2  (4.3- 13a) 

where 
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(4.3- 13b) 


c  =  inf  gfR^Cjco)] 

<?(A)  is  the  minimum  singular  value  of  A. 

Theorem  4.3-6  is  very  useful  in  speeding  up  the  convergence  rate.  Notice  that 
Theorem  4.3-4  and  Corollary  4.3-5  are  so  alike  that  the  only  difference  is  the  direction  of 
the  inequality  arrows.  Corollary  4.3-5  can  be  obtained  by  reversing  the  inequality  arrows  in 
Theorem  4.3-4.  This  inequality-arrow-reversing  property  motivated  us  that  Theorem  4.3-6 
may  have  its  counterpart  as  follows. 

Theorem  4.3-7: 

If  p(y)  <  1,  then 

y0  *  -J Ji2(y)  y2  +  [  1  -  M-2(Y)  ]  c2  (4.3-14) 

where  c  is  given  by  (4.3- 13b). 

Proof :  See  Appendix  A. 

The  following  corollary  plays  an  important  role  in  our  algorithm.  This  corollary  can 
be  easily  derived  by  using  the  properties  in  Theorem  4.3-1. 

Corollary  4.3-8: 

Given  y{,  yv  M-(YX)-  and  p(y2)  and  yx  <  y2-  The  optimal  norm  yo  lies  between  yl 
and  y2,  i.e.,  Yj  <  yo  <  yr  Suppose  y5  is  a  point  such  that  yx  <  y5  <  yv  and  the 
value  p(y5)  is  known. 

i)  If  p(y5)  >  1,  then 

Y0  >  Y5  +  [H(Y5)  -  i][y5  -  yx\  /  tiKYj)  -  n(y5)]  (4.3-15) 

and 

Y0  <  Y5  +  tH(Y5)  -  1][Y5  -  Y21  /  tM-(Y2)  *  M-(Y5)1  (4.3-16) 

ii)  If  p(y5)  <  1,  then 

Y0  <  Y5  +  ^-IMYs-YJ/OKYj)-^)]  (4.3-17) 
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and 


Y0  >  Y5  +  &KY5)  -  ^fYj  -  Y21  /  DKY2)  -  JKYj)] 


(4.3-18) 


Proof  of  Corollary  4.3-8:  See  Appendix  A. 

4.3.3  Fast  y-Iteration  Algorithm 

We  will  start  with  the  initial  lower  and  upper  bounds  of  Y-  The  initial  lower  and 
upper  bound  are 

YL  =  max  {  Yl1.Yl3  }  (4.3-19) 

and 

Yu  =  tnin  {  Ym,  YU2  )  (4.3-20) 

where  Yuj.Yu2.Ylj.  and  yL3  are  given  by  (4.3-5a),  (4.3-5b),  (4.3-6a),  and  (4.3-6c) 
respectively. 

Before  entering  the  iteration  loop,  some  initialization  should  be  done.  Let 

Yj  =  Yl  and  y2  =  yv  (4.3-21) 

and  evaluate  pty)  and  p(Y2)-  Since  Yj  <  Yc  <  Y2»  we  have 

p(Yj)  >  1  and  p(y2)  <  1  (4.3-22) 

From  Theorem  4.3-4,  Theorem  4.3-6,  Corollary  4.3-5,  and  Theorem  4.3-7,  we 
have  the  following: 


VI 

y^2(Yj)Yj+[i 

H2(Y,)]  c2 

•'=  <*U1 

(4.3-23) 

Yo  <  , 

-  M.2(Y2)3  b2 

:=  au2 

(4.3-24) 

Yo  >, 

■  M'2(Y1)]  b2 

:=  «L2 

(4.3-25) 

and 
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(4.3-26) 


y0  *  yJvHfj  f2  +  [l-  M-2(Y2)]  c2  :=  aLi 
where  b  =  II  R2i  It  and  c  =  inf  a[R,.(jco)]. 

00  CO  - 

It  is  easy  to  o&e  that  the  intersection  of  the  horizontal  line  y  =  1  and  the  straight  line 
connecting  the  two  points  (Yrli(y1))  and  (Y2.ti(Y2»  is  at  the  right  of  the  point  (yo,1).  That 

is. 


Y0  <  Y2  +  lti(Y2) '  1][Y2-  Yj  /  [M-(Yj)  -  ti(Y2)l  (4.3-27) 

Let 

Y3  =  max  (ctn.c^)  (4.3-28) 

Y4  =  min  {  am,  ocU2,  am  )  (4.3-29) 

where  0^,  aLV  0^,  a^,  and  are  given  by  (4.3-23)  -  (4.3-27).  It  is  easy  to  see  that 

Yi  <  Y3  ^  Y0  ^  Y4  <  Y2  (4.3-30) 

After  the  initialization  of  Yi,  ti(Yi).  Y2.  ti(Y2)>  Y3>  and  Y4,  we  are  ready  to  enter  into 
the  iteration  loop.  The  first  step  of  the  iteration  loop  is  to  let 

Y5  =  (Y3  +  Y4)/2  (4.3-31) 

and  compute  p(Y5).  Define 

Pi  :=  7li2(Y5)^  +  [l-li2(Y5)]c2  (4.3-32) 

P2  :=  M-2(Y5)  Y^  +  t1  -  M-2(Y5)]  fc2  (4.3-33) 

P3  :=  Y5  +  tti(Y5)  -  IKY5 '  Yj  /  WV  -  ti(Y5)l  (4.3-34) 

P4  :=  Y5  +  tH(Y5)  -  1][Y5  -  Y2]  /  fc(Y2)  *  H(Y5)1  (4.3  35) 

where  b  =  II  R21 II  and  c  =  inf  a[R,.(ja))]. 

Oft 
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Now  we  compare  |x(y5)  with  1.  If  ji(y5)  >  1,  then  from  Theorem  4.3-6,  Corollary 
4.3-5,  and  Corollary  4.3-8,  we  have 

Yo  ^  Pi.  Y0  >  P2.  Y0  >  Ps.  Y0  <  P4  (4.3-36) 

Thus,  y3,  y4,  yt,  and  pty)  are  updated  as  follows 

Y3  =  max  {  p2,  P3  },  Y4  -  min  {  fy,  p4,  y4)  , 

Yj  =  Y5 .  M-fy)  =  H(Y5)  (4.3-37) 

where  pj,  P2,  P3,  and  P4  are  given  by  (4.3-32)  -  (4.3-35). 

If  H(Y5)  <  1.  then  from  Theorem  4.3-4,  Theorem  4.3-7,  and  Corollary  4.3-8,  we 

have 

Yo  <  P2 .  Yo  ^  Pi .  Yo  <  Ps .  Yo  >  P4  (4.3-38) 

Thus,  y3,  y4,  y2,  and  p.(y2)  are  updated  as  follows 

y3  =  max  {  pr  p4,y3  },  Y4  =  min  {  p2,  P3}  , 

Y2  =  Y5 .  H(y2)  -  n(y5)  (4.3-39) 

where  Pr  P2,  P3,  and  P4  are  given  by  (4.3-32)  -  (4.3-35). 

The  iteration  loop  is  repeated  until  the  difference  y4  -  y3  is  small  enough.  The 
accuracy  of  yo,  i.e.,  ( y4  -  y3  )  /  2,  can  be  arbitrarily  small  by  simply  increasing  the  number 
of  iterations. 

The  fast  iterative  algorithm  for  the  computation  of  yo  in  (4.3-3)  is  summarized  as 
follows. 

Fast  y-iteration  Algorithm: 

(1 )  Compute  the  initial  lower  and  upper  bounds 

Yj  =  max  {  yu,  yu  ),  y2  =  min  {  yuit  yu2  ) 
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where  Yur  yU2,  yl1,  and  yu  are  given  by  (4.3-5a),  (4.3-5b),  (4.3-6a),  and  (4.3- 

6c)  respectively. 

(2)  Evaluate  jx(Yj)  and  |i(y2). 

(3)  If  JJ-CYj)  or  (x(Y2)  equals  to  1,  then  Y0  =  I^(Y1)  or  Y0  =  H(Y2)  and  go  to  (11). 

(4)  Compute  the  bounds  y3  and  y4. 

Y3  =  max  {  c^,  au  },  Y4  =  min  {  am,  am,  aU3  } 
where  c^,  a^,  auv  a^,  and  aU3  are  given  by  (4.3-23)  -  (4.3-27). 

(5)  If  I Y4  *  Y3 1  <  £.  then  Y0  =  ( Y3  +  Y4 )  /  2  and  go  t  o  (11).  (e  is  the  acceptable 
accuracy.) 

(6)  Let  Y5  =  (Y3  +  Y4)/2  and  evaluate  |1(y5). 

(7)  If  ji(y5)  =  1,  then  yo  =  Y5  and  go  to  (11). 

(8)  If  }!(y5)  <  1,  then  go  to  (10). 

(9)  If  ji(y5)  >  1,  then  update  yy  y4,  Yt  and  jity)  as 

Y3  =  max  {  P2,  P3  },  Y4  =  tnin  {  pr  p4,  y4} 

Yx  -  Y5  .  li(Yj)  =  H(YS) 

and  then  go  to  (5),  where  pr  p2,  p3,  and  p4  are  given  by  (4.3-32)  -  (4.3-35). 

(10)  Update  yy  y4.  Y2  and  p(y2)  as 

Y3  =  max  {  pj,  p4,  y3  },  Y4  =  min  {  p2,  P3  } 

Y2  =  Y5  .  \i(y2)  =  4(Y5) 

and  then  go  to  (5),  where  pj,  P2,  p3,  and  p4  are  given  by  (4.3-32)  -  (4.3-35). 

(11)  Stop. 
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4.4  Using  DGKF's  Solution  To  Solve  A  Four-Block  H~  Optimization 
Problem 


In  [29],  DGKF  assume  the  realization  of  G(s)  is  given  by 


’  A 

Bi 

V 

II 

C/3 

''W' 

O 

Ci 

0 

D12 

C2 

D21 

0 

with  the  following  assumptions. 

(i)  (A,  B^  is  stabilizable  and  (Clt  A)  is  detectable. 

(ii)  (A,  B2)  is  stabilizable  and  (C2,  A)  is  detectable. 

(iii)  d[2 [ Ci  D12]  =  [ 0  I]. 


■  B1 " 

II 

Hd 

O' 

_  D21  _ 

21 

.  I . 

Define  two  Hamiltonian  matrices  as  follows. 


Ho„  := 


T 

-C  C 
uri 


Y\BTrB2BT2 


-A* 


and 


:= 


AT  r2cfa-CjC2 

-B|B]  -A 


(4.4-1) 


(4.4-2a) 

(4.4-2b) 

(4.4-2c) 

(4.4-2d) 


(4.4-3a) 


(4.4-3b) 


Then  the  following  theorems  will  characterize  the  set  of  suboptimal  stabilizing  controllers 
such  that  IIOII..  <  y  where  C>  is  the  closed-loop  transfer  matrix  from  v  to  z. 
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Theorem  4.4-1:  [29] 

There  exists  a  stabilizing  controller  such  that  It  IL  <  Y  if  and  only  if  the  following 
three  conditions  hold. 

(i)  HL.  6  dom(Ric)  and  :=  Ric(HJ  £  0.  (4.4-4a) 

(ii)  €  dom(Ric)  and  Y„  :=  Ric(J,J  £  0.  (4.4-4b) 

(iii)  p(X00YJ<y2.  (4.4-4c) 

Moreover,  when  these  conditions  hold,  one  such  controller  is 

(4.4-5a) 

where 

K  =  A  +  Y2B,bX  +  +  Z_LC2  (4.4-5b) 

F„  :=  -B2X_ ,  U.:=-YMC2,  Z»  :=  (I -Y2YJ^yl .  (4.4-5c) 


In  the  above  theorem,  condition  (i)  means  that  there  exists  a  nonnegative  definite 
solution  X„  to  the  Riccati  equation  corresponding  to  the  Hamiltonian  H^.  Condition  (iii) 
means  that  the  spectral  radius  of  X„YM  is  less  than  t2. 

Theorem  4.4-2:  [29] 

If  conditions  (i)  to  (iii)  in  Theorem  4.4-1  are  satisfied,  the  set  of  all  stabilizing  controllers 
such  that  II  <X>  IL<  y  equals  the  set  of  all  transfer  matrices  from  y  to  u  in 


where  Q  e  RH~,  IIQIL<y. 


(4.4-6) 


The  above  two  theorems  show  an  easy  state-space  approach  in  constructing  a 
stabilizing  subcptimal  controller  such  that  IIOIL  <  Y-  The  order  of  the  suboptimal  controller 
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can  be  the  same  as  that  of  the  plant  G(s).  The  major  computation  involved  is  the  solution  of 
two  Riccati  equations  which  are  easy  to  solve  if  solution  exist. 

The  DGKF  approach  is  a  great  break-through  in  the  solution  of  H°°  optimization 
problem.  However,  if  we  try  to  reduce  the  value  of  y  to  the  minimum  we  will  encounter 
numerical  difficulty  before  y  reaches  the  minimum.  The  optimum  occurs  when  y2  equals  to 
pCX^Y^).  When  y2  is  close  to  p(X00Y00),  the  matrix  (I  -  Y2Y00X00)  becomes  ill-conditioned 
and  the  inversion  =  (I  -  Y2YaoX^)‘ 1  will  not  be  numerically  reliable.  In  (4.4-5a),  both  A 
and  B  matrices  of  Ksub(s),  i.e.,  and  -ZJ^^,  rely  on  the  computation  of  Zm  and  therefore 
on  the  the  inversion  of  (I  -  Y2YJX^). 

In  [38],  we  attempted  to  use  DGKF  approach  to  design  an  optimal  controller  for  a 
four-block  H“  optimization  problem.  We  started  from  a  large  Y  =  yu  which  guarantees 
solution  existence  of  Riccati  equations  and  pCX^Y,.)  <  y2.  Then  y  is  updated  as  y  <— 
V  P(X«,YJ ,  the  new  Riccati  equation  solutions  X,*,  and  YM  will  make  pfX^Y.J  >  y2.  Let  y^ 
<-  y  and  y  «-  V  P(X«Y.J,  the  new  Riccati  equation  solutions  X„  and  Y„  will  make 
pfX^Y..)  <  y2.  Let  y0  <—  y  and  y  «-  •/  p(X^Y„).  The  iterative  process  is  repeated  until  the 
gap  Yu  -  yL  is  small  enough.  The  convergence  rate  is  fast.  For  many  examples,  it  only  takes 
8  to  10  iterations  to  reduce  the  gap  from  100  to  10'10. 

As  y  is  close  to  the  optimum,  i.e.,  the  gap  Yu  -  Ji  is  approaching  to  zero,  the  matrix 
(I  -  Y2YwXoc)  will  be  nearly  singular  and  then  the  elements  of  and  A^  will  approach  to 
infinity.  It  looks  like  that  there  is  no  way  to  design  a  controller  such  that  IIOIL  is  close  to  the 
optimum  since  the  elements  of  A  and  B  matrices  of  the  controller  will  all  approach  to 
infinity.  Nevertheless,  the  gain  of  the  controller  is  not  really  infinite.  If  we  do  partial 
fraction  expansion  for  the  controller  we  will  find  that  only  some  of  the  controller  poles  and 
their  corresponding  gains  will  approach  to  infinity  and  the  rest  will  remain  finite  when  y  is 
close  to  the  optimum.  These  terms  with  infinite  poles  and  gains  can  be  approximated  by 
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finite  directly  feedthrough  terms.  In  other  words,  the  practical  optimal  controllers  will  be 
proper  with  direct  feedthrough  terms  instead  of  strictly  proper. 

The  numerical  difficulty  is  mainly  caused  by  the  restriction  of  the  controller  to  be 
strictly  proper.  Most  of  optimal  H~  controllers  are  proper  with  direct  feedthrough.  If  a 
proper  with  direct  feedthrough  optimal  controller  is  forced  to  be  represented  in  a  strictly 
proper  form,  some  poles  and  their  corresponding  gains  will  be  forced  to  be  infinity. 

The  following  is  a  simple  four-block  H”  optimization  problem  which  is  used  to 
illustrate  our  numerical  experience  in  using  the  DGKF's  approach  to  design  an  H“  optimal 
controller.  A  realization  of  the  plant  G(s)  is  given  by 


'  A 

B1 

B.  " 

'-1 

0 

1 

0 

o’ 

2 

0 

2 

0 

0 

1 

G(s)  = 

D11 

°12 

= 

1 

1 

0 

0 

0 

(4.4-7) 

0 

0 

JL 

0 

1 

C2 

D21 

D22_ 

_  1 

1 

0 

i 

0 

The  assumptions  in  (4.4-2)  are  all  satisfied  except  assumption  (i)  which  we  believe  is  too 
restrictive.  However,  Theorems  4.4-1  and  4.4-2  are  still  applicable  to  this  example  even 
(A,  B,)  is  not  stabilizable. 

We  started  from  y  =  100.  Conditions  (i)  to  (iii)  are  all  satisfied  and  7  pQOO 
=  4.647641 1024  which  is  less  than  y  =  100.  Update  y  to  be  4.647641 1024  and  compute 
the  new  Riccati  equation  solutions  X*.  and  YM  and  7  P(X~Y00)  =  4.7375420092  which  is 
greater  than  y  =  4.6476411024  and  condition  (iii)  is  not  satisfied.  The  new  value  of 
7  P(X<~YJ  ,  i.e.,  4.7375420092  is  used  to  update  y.  Now,  just  two  iterations  after  the 
initial  guess  y  =  100,  we  find  that  the  optimal  IIOIL,  y0,  is  between  4.6476411024  and 
4.7375420092.  After  six  more  iterations,  we  have 

4.7341604761  <  y0  <  4.7341604768  <4-4'8) 


52 


With  y  =  4.7341604768  which  is  very  close  to  the  optimum,  from  (4.4-5)  we  have  a  state- 
space  realization  of  the  controller  as  follows: 


with 


K(s)  = 


2.77 85374772e+09  2.7785374782e+09 

-2.4722482140e+10  -2.4722482142e+10 


(4.4-9a) 


(4.4-9b) 


-2.77 85374782e+09 
2.4722482140e+10 


(4.4-9c) 


^  =  [-3.10880347:86-01  -4.2370320107e+00  ]  (4.4-9d) 


The  eigenvalues  of  AK  are  -8.7542343 140e-01  and  -2.1943944664e+10.  Notice  that  one 
of  the  controller  poles  is  pretty  far  away.  When  y  is  nearly  optimal,  this  pole  will  move 
further  to  the  infinity.  Do  a  partial  fraction  expansion  for  the  above  second  order  strictly 
proper  controller  K(s),  we  have 


K(s)  = 


-8.9043458823e^  jO 
s+8.7542343 140e-01 


+ 


1.0388615552e+ll 

s+2.1943944664e+10 


(4.4-10) 


The  second  term  of  the  above  expression  is  a  wide  band  low-pass  filter  which  can  be 
approximated  by  a  direct  feedthrough  term;  therefore,  K(s)  can  be  approximated  by  the 
following: 


K  ,  »  _  -8.9Q43458823e+OQ 
~  s+8.7542343 140e-01 


+  4.7341604762 


(4.4-11) 


With  the  controller  K^s),  the  closed-loop  system  is  internally  stable  and  the  norm  of  the 
closed- loop  system  from  v  to  z,  NOIL  is  4.7341604762  which  is  very  close  to  the  optimum 
up  to  the  accuracy  of  10'10.  This  nearly  optimal  controller  is  proper  with  direct  feedthrough 
and  is  first  order  which  is  one  order  less  than  that  of  the  plant. 
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SECTION  5 


COMPUTATION  OF  THE  REAL  STRUCTURED  SINGULAR  VALUE 

Stability  robustness  is  an  important  issue  in  the  analysis  and  design  of  control 
systems.  Currendy,  there  are  two  major  approaches  to  stability  robustness  analysis.  One  is 
the  structured-singular-value  (SSV)  [39-45]  or  the  multivariable-stability-margin  (MSM) 
[46-50]  approach  and  the  other  is  the  perturbed-characteristic-polynomials  approach  [51- 
56].  Several  significant  progresses  have  been  made  in  both  approaches. 

For  a  large  class  of  linear  time-invariant  systems  with  real  parametric  perturbations, 
the  coefficient  vector  of  the  characteristic  polynomial  is  a  multilinear  function  of  the  real 
parameter  vector.  Based  on  this  multilinear  mapping  together  with  the  recent  results  by  De 
Gaston  and  Safonov  [48],  Sideris  and  Pena  [50],  Bartlett,  Hollot,  and  Lin  [52],  and 
Bouguerra,  Chang,  Yeh,  and  Banda  [57],  an  algorithm  for  computing  the  real  structured 
singular  value  is  proposed.  The  algorithm  requires  neither  frequency  search  nor  Routh's 
array  symbolic  manipulations  and  allows  the  dependency  among  the  elements  of  the 
parameter  vector.  Moreover,  the  number  of  the  independent  parameters  in  the  parameter 
vector  is  not  limited  to  three,  as  is  required  by  many  existing  structured  singular  value 
computation  algorithms. 

In  Sec.5.1,  the  definition  of  the  real  structured  singular  value  is  reviewed,  and  the 
basic  concepts  of  the  proposed  algorithm  are  briefly  described.  In  section  5.2,  the  mapping 
properties  from  the  parameter  space  to  the  coefficient  space  and  the  theories  on  which  the 
iterative  algorithm  is  developed  will  be  demonstrated.  The  detailed  iterative  algorithm  for 
computing  the  real  structured  singular  value  is  summarized  in  Sec.5.3.  Sec.5.4  is  a 
summary  of  the  fast  segment  stability  checking  algorithm  which  is  used  repeatedly  in  the 
proposed  iterative  algorithm  in  computing  the  real  SSV. 
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5.1  Definition  of  the  Real  Structured  Singular  Value  and  Preliminaries 

According  to  Doyle  and  Safonov  [39-40,46-47],  all  the  plant  uncertainties, 
structured  or  unstructured,  unmodeled  dynamics  or  parametric  perturbations,  can  be 
described  by  the  block  diagram  shown  in  Fig.5.1-1.  In  Fig.5.1-1,  A(s)  =  block  diag 
{ A^s),  A2(s),  ...,  Am(s)}  and  M(s)  is  the  nominal  system  which  includes  the  nominal 

plant  and  the  stabilizing  controller.  Doyle  [39-40]  and  Safonov  [46-47]  defined  the 
structured  singular  value  (SSV)  and  the  multivariable  stability  margin  (MSM)  respectively 
based  on  the  above  perturbation  structure  and  used  them  as  analysis  tools  for  robust 
stability.  The  SSV  and  the  MSM  are  nonconservative  scalar  stability-margin  measures  for 
multivariable  systems. 


Fig.5.1-1  Standard  structure  for  a  perturbed  closed-loop  system. 

Algorithms  [39-45]  to  compute  the  SSV  are  available  only  for  those  cases  where  the 
number  of  perturbation  blocks  are  less  than  or  equal  to  three.  The  computational  problem 
for  the  cases  with  more  than  three  perturbation  blocks  is  still  an  unsolved  problem. 

One  important  special  case  of  plant  uncertainties  is  the  real  parametric  perturbation. 
In  this  case  the  perturbation  matrix  A(s)  is  a  real  diagonal  matrix.  The  SSV  defined  for  this 
case  is  called  the  real  SSV.  Algorithms  of  computing  the  real  SSV  [44-45]  are  also  only 
available  for  those  cases  where  the  number  of  independent  perturbation  parameters  are  less 
than  or  equal  to  three.  The  MSM  defined  for  this  case  is  called  the  real  MSM.  An  iterative 
algorithm  of  computing  the  real  MSM  for  real  diagonal  A  was  developed  by  De  Gaston  and 
Safonov  [48]  and  generalized  by  Pena  and  Sideris  [49].  There  is  no  limitation  on  the 
number  of  perturbation  parameters.  However,  this  iterative  algorithm  is  complicated  since 


an  extensive  frequency  search  is  required.  In  [50],  Sideris  and  Pena  eliminated  the  need  of 
frequency  search  by  requiring  the  first  column  of  the  Routh's  array  to  be  positive.  This 
approach  requires  symbolic  manipulations.  Besides,  the  first-column  elements  of  the 
Routh's  array  are  not  multilinear  functions  of  the  parameter  vector  even  though  the  original 
characteristic  coefficients  are.  For  those  elements  to  be  multilinear  functions  of  the 
parameters,  more  dependent  parameters  need  to  be  created  which  will  cause  unnecessary 
complexity. 

In  the  following  we  assume  that  the  perturbation  matrix  A  in  Fig.5.1-1  is  real 
diagonal,  i.e.,  A  =  diag  {  Sj,  §2,  ...,  Bm  }  and  the  nominal  system  M(s)  is  a  rational  matrix 
with  real  coefficients.  If  the  parameters  vary  independently  and  -1£8.  £  1,  i  =  l,2,...,m, 

the  parameter  perturbation  domain  JS  can  be  described  as  a  hyper-cube  JS)  with  2m  vertices 
(±1, ...,  ±1)  in  the  m-dimensional  real  space.  In  general,  the  perturbation  matrix  A  can  be 
written  as  A  =  diag  {8jlml,  82Im2, ...,  8^^}  where  1^  is  the  identity  matrix  of  order  mi 
and  ml+m2+  ...+mr  =  m.  That  is,  5=BL=  •••  =8  =8,,  8  ,  =8  ,  „=  =8  ,  „=8„, 

’12  ml  r  ml+1  ml +2  ml+m2  2’ 

...,  etc.  In  this  case,  the  parameter  perturbation  domain  JS  is  an  r-dimensioral 

hyperrectangle  inside  the  m-dimensional  hypercube  JS.  The  system  is  said  to  be  robustly 
stable  in  JS  if  and  only  if  it  is  stable  for  every  parameter  vector  8  =  [  8j  82  ...  8m  ]T  in 

JS.  Throughout  the  report,  we  may  use  "JS  is  stable"  to  replace  "the  system  is  robustly 
stable  in  JS"  whenever  there  is  no  confusion. 

The  real  multivariable  stability  margin  (real  MSM)  kj^  is  defined  as  the  largest  real 

constant  k  such  that  the  closed-loop  system  remains  robustly  stable  in  kJS  where  kJS  is  the 
enlarged  (or  shrunk)  parameter  perturbation  domain  of  J9,  i.e., 

kJS  :=  {  8:  8  =  [Sj.-.Sj,  82..,82, ....  8r,..,8r]  e  Rm 

and  ISJSk,  i=l,2,...,r  }  (5.1-1) 

The  enlarged  (or  shrunk)  hypercube  of  JS,  kJS  is 

kJS  :=  {  8:  Be  Rm  and  I  8.1  £  k,  i=l,2,...,m  }.  (5.1-2) 
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Recall  that  the  real  structured  singular  value  (real  SS  V)  p  is  defined  as 

ji  :=  [  min  {  k  I  det  [I+M(jco)A]  =  0  for  some  co  and  As  X  (k)  }  ]  1  (5.1-3) 

where 

X  (k)  =  {  A  I  diag  {8^,  8^ . 8^}  with  I  8j  I  <;  k,  for  all  i }  (5.1-4) 

That  is,  the  real  SSV  p  is  the  reciprocal  of  the  smallest  k  such  that  the  system  is  unstable  in 
kJ9.  It  is  easy  to  see  that  the  relation  between  the  real  SSV  p  and  the  real  MSM  kM  is 

p  =  1  /  kM  (5.1-5) 

Several  significant  results  [51-56]  have  been  obtained  in  the  perturbed- 
characteristic-polynomials  approach.  Probably  the  most  famous  are  the  Kharitonov's 
theorems  [51]  which  apply  to  the  special  case  with  a  hyper-rectangular  perturbed  region  in 
the  coefficient  space.  In  this  special  case,  the  coefficients  of  the  characteristic  polynomial 
vary  independently  and  the  robust  stability  of  the  system  can  be  easily  determined  by  four 
bounding  characteristic  polynomials.  Unfortunately,  the  Kharitonov's  theorems  cannot  be 
applied  to  our  problem  since  the  coefficient  variations  of  the  characteristic  polynomial  are 
not  independent. 

Bartlett,  Hollot,  and  Lin  [52]  developed  an  important  theorem  which  is  applicable  to 
the  case  when  the  coefficients  of  characteristic  polynomial  are  linearly  dependent  The 
theorem  is  now  well  known  as  the  Edge  Theorem:  For  the  set  of  characteristic  polynomials 
inside  a  polytope  IP  in  the  coefficient  space,  every  polynomial  in  IP  is  stable  if  and  only  if 
all  the  exposed  edges  of  IP  are  stable.  This  simplifies  the  stability  checking  tremendously 
since  checking  the  stability  of  exposed  edges  is  much  simpler  than  checking  that  of  the  full 
IP.  The  exposed-edge  stability  checking  is  done  by  sweeping  t  from  0  to  1  such  that 

ta‘  +  (l-t)aj  (5.1-6) 

are  all  stable  for  all  vertices  a1  and  a*  of  F. 
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bialas  [54]  and  Fu  and  Barmish  [53]  reduced  the  checking  of  the  exposed-edge 
sweep  stability  to  a  one-shot  test.  They  showed  that  ta‘+  (1-t)  a-"  is  stable  for  all  t  e 
[0,1]  if  and  only  if  the  real  eigenvalues  of  -H.H.'1  are  all  negative  where  a1  and  a-*  are 
assumed  to  be  stable  and  H.  and  H  are  the  Hurwitz  matrices  for  a'  and  aJ  respectively. 

Recently,  a  fast  algorithm  based  on  Chapellat  and  Bhattacharyya's  Segment  Lemma  [58] 
was  proposed  by  Bouguerra,  Chang,  Yeh,  and  Banda  [57]  for  checking  the  stability  of  the 
exposed  edges.  The  computation  in  the  algorithm  mainly  depends  on  the  number  of  vertices 
instead  of  the  edges  and  therefore  reduces  the  computation  burden  due  to  the  "combinatoric 
explosion." 

There  are  no  such  celebrated  properties  in  the  parameter  space  as  those  in  the 
coefficient  space  discovered  by  Kharitonov  [51],  Bartlett,  Hollot,  and  Lin  [52],  The 
closed-loop  system  may  be  unstable  inside  J§  although  it  is  stable  at  all  the  edges  of  the 
hypercube  JS>  [59].  So  far,  there  is  no  easy  way  of  checking  robust  stability  in  the 
parameter  space. 

For  each  parameter  vector  8  in  the  parameter  perturbation  domain  J9,  there  is  a 
corresponding  characteristic  polynomial,  i.e.,  a  coefficient  vector  a  in  the  coefficient 
space.  Let  Jl(J9)  be  the  image  of  J9  in  the  coefficient  space.  The  closed-loop  system  is 
robustly  stable  in  JD  if  and  only  if  every  characteristic  polynomial  in  d(J$>)  is  stable. 
Although  several  significant  results  for  robust  stability  have  been  obtained  in  the  coefficient 
space,  there  is  no  efficient  way  to  check  robust  stability  for  $(«©)  since  d(J3)  usually  is 
neither  a  Kharitonov's  hyper-rectangle  [51]  nor  a  poly  tope  considered  by  Bartlett,  Hollot, 
and  Lin  [52]. 

Define  the  polytope  P(J§)  as  the  convex  hull  of  the  2m  image  points  in  the  (n+1)- 
dimensional  coefficient  space  mapped  from  the  vertices  of  J9  where  n  is  the  degree  of  the 
characteristic  polynomial.  If  the  mapping  is  multilinear,  then  the  image  of  the  edges  of  the 
hypercube  J9  will  be  the  straight  line  segments  connecting  the  corresponding  mapped 
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vertices.  The  image  of  J0,  d(,E>),  is  a  subset  of  d(J9)  and  therefore  is  a  subset  of  the 
polytope  F(J9).  Under  the  condition  of  the  multilinear  mapping,  the  stability  of  the  edges 
of  1P(JS)  will  guarantee  the  robust  stability  in  JS>.  The  multilinear  mapping  can  be  easily 
achieved  by  assuming  that  the  nominal  system  M(s)  in  Fig.5.1-1  is  a  rational  matrix  with 
real  coefficients. 

Now,  we  have  an  easy  way  to  check  the  sufficient  condition  for  the  robust  stability 
in  J 0  by  using  its  corresponding  polytope  1P(JS).  The  sufficient  condition  is  still  not 
enough  to  determine  the  real  MSM  kM.  Any  k  such  that  the  polytope  IP(kjjj)  remains 
stable,  say  kL,  can  be  served  as  a  lower  bound  for  kM.  However,  there  may  exist  some  k  > 
kL  such  that  kJS)  is  stable  although  the  corresponding  polytope  !P(kji9)  is  unstable. 

Any  k  which  cause  instability  of  kJS  or  d(kJ0)  qualifies  as  an  upper  bound  for  kM. 

To  facilitate  the  description  of  the  relations  between  the  parameter  space  and  the  coefficient 
space,  let  the  edges  (vertices,  resp.)  of  ^(kJS)  which  are  mapped  from  the  edges  (vertices, 
resp.)  that  are  parallel  to  an  axis  of  coordinates  of  kJS)  be  called  the  crucial  edges  (vertices, 
resp.)  and  those  which  are  not  crucial  be  called  noncrucial  edges  (vertices,  resp.).  The 
noncrucial  edges  include  two  kinds  of  edges:  supplemental  edges  and  fictitious  edges.  The 
supplemental  edges  are  the  image  of  the  edges  of  kj0  which  are  not  in  kJS.  The  fictitious 
edges  are  the  edges  of  1P(kJS>)  which  are  not  mapped  from  the  edges  of  kJ9.  The  crucial 
edges  are  all  in  d(kJ9).  Thus,  some  k  which  causes  instability  at  the  crucial  edges  of 
IP  (k  J9),  say  ky,  may  be  used  as  an  upper  bound  for  kM.  If  the  lower  and  the  upper  bounds 

coincide  or  are  close  enough,  we  have  the  real  MSM  kM  and  the  real  SSV  (i.  The  objective 

of  the  iterative  algorithm  is  to  reduce  the  gap  of  the  lower  and  the  upper  bounds.  When  the 
gap  is  smaller  than  the  desired  accuracy  e,  i.e.,  Iky  -  kj  $  e,  we  have  the  real  MSM  kM=  kL 
and  the  real  SSV  p  =  l/kM. 
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5.2  Multilinear  Mapping  and  the  Polytopic  Polynomials 

The  robust  stability  of  the  perturbed  system  in  Fig.5.1-1  is  determined  by 
det[I+M(s)A].  There  are  two  lemmas  available  for  checking  the  robust  stability.  They  are 
listed  as  follows. 

Lemma  5.2-1: 

The  closed-loop  system  is  robustly  stable  in  29  if  and  only  if  M(s)  is  asymptotically 
stable  and 

det  [  I  +  M(jco)A  ]  *  0  for  all  co  and  all  8  in  29.  (5.2-1) 

Lemma  5.2-2: 

The  closed-loop  system  is  robustly  stable  in  29  if  and  only  if  all  the  zeros  of 

det  [  I  +  M(s)A  ]  for  all  S  in  29  (5.2-2) 

are  in  the  open  left  half  of  s-plane. 

The  existing  computational  algorithms  for  the  SSV  and  the  MSM  are  all  developed 
based  on  Lemma  5.2-1  while  the  approach  to  be  presented  in  Sec.5.3  is  based  on  Lemma 
5.2-2. 

Recall  that  A  =  diag  {5^,  S2Im2, ....  5^}  =  diag  {8r  S2 . 8J  where  1^  is 

the  identity  matrix  of  order  mi  and  ml+m2+  ...+mr  =  m.  That  is,  5j=82=  •••  =5ml=8r 

^mi+i=^mi+2=  =^mi+m2=^2’  **”  etc-  ^he  par*111®1®1,8  82,  ...»  8f  are  assumed 

independent  and  -1  ^  8.  ^  1,  i  =  1,2, ...,r.  The  parameter  perturbation  domain  29  is  an  r- 

dimensional  hyperrectangle  inside  the  m-dimensional  hypercube  J9  with  2m  vertices  (±1, 
...,  ±1)  in  the  m-dimensional  real  space. 

First  of  all,  we  will  establish  the  relationship  between  the  coefficients  {a0,  ar  o^, 
...,  an)  of  the  characteristic  polynomial 

a0sn  +  OjS11'1  +  a2sn'2  + . +  an  (5.2-3) 
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and  the  perturbation  parameters.  It  is  interesting  to  note  that  if  M(s)  is  a  rational  matrix  with 

real  coefficients  then  the  coefficient  vector  a  =  [  aQ  aj  a2  ...  an  ]T  is  a  multilinear 

function  of  the  parameter  vector  8  =  [  8}  S2  ...  8m  ]T.  That  is,  if  we  only  allow  one  of  the 

B.'s,  say  BL,  to  vary  and  keep  the  rest  m-1  of  the  5/s  constant  then  a  is  a  linear  function  of 

8.. 

j 

Theorem  5.2-3: 

Consider  the  perturbed  system  in  Fig.5.1-1  where  A  =  diag  {  8j,  82,  ...,  8m  }  and  M(s) 
is  rational  with  real  coefficients.  The  coefficients  a/s,  i=0,l,2,...,n,  of  the 
characteristic  polynomial  of  the  perturbed  system  are  multilinear  functions  of  8.'  s . 

Proo'.  See  Appendix  A. 

Remark  5.2-4: 

The  coefficients  a.'s,  i=0,l,2,...,n,  of  the  characteristic  polynomial  of  the  perturbed 
system  in  general  are  not  multilinear  functions  of  8/s  unless  JS>  =  3. 

The  following  lemma  is  a  direct  consequence  of  the  multilinear  mapping  between  8 

and  a. 

Lemma  5.2-5: 

In  the  parameter  space,  draw  a  line  segment  with  ending  points  Ej  and  E2  and  the  line  is 

parallel  to  an  axis  of  coordinates.  The  images  of  these  two  ending  points  are  denoted  by 
d(Ej)  and  d(E2)  respectively.  The  image  of  the  line  segment  in  the  coefficient  space  is  a 
straight  line  segment  which  connects  d(Ej)  and  d(E2)  if  the  mapping  from  the  parameter 

space  to  the  coefficient  space  is  multilinear. 

The  parameter  perturbation  domain  J9  is  an  r-dimensional  hyperplane  inside  the  tri¬ 
dimensional  hypercube  3  with  2m  vertices  whose  edges  are  parallel  to  the  axes  of 
coordinates.  d(J9)  and  t)(JS)  are  the  images  of  JS  and  3  respectively  in  the  coefficient 
space.  The  polytope  P  (3)  is  the  convex  hull  of  the  2m  image  points  in  the  (n+1)- 
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dimensional  coefficient  space  mapped  from  the  vertices  of  3.  In  the  following,  we  will 
show  that  if  the  mapping  is  multilinear,  the  image  of  3,  i.e.,  $(JS),  is  a  subset  of  the 
polytope  1P(JS>). 

Theorem  5.2-6: 

d(J)  <=  (5.2-4) 

Proof:  See  Appendix  A. 

Theorem  5.2-7: 

If  c  i>2,  then 

d(J®j)  <=  $<JS2)  (5.2-5) 

and 

POSj)  <=  F(J92).  (5.2-6) 

Proof:  See  Appendix  A. 

Theorem  5.2-8: 

The  hypercube  or  hyperrectangle  3  is  cut  into  two  equal  subdomains  3 1  and  32  by  a 
hyperplane  which  is  orthogonal  to  an  axis  of  coordinates.  Then 

4(JS)  <=  F(^j)uF(J®2)  c  Jp(3)  (5.2-7) 

Proof:  See  Appendix  A. 

In  the  following,  a  simple  case  with  2-dimensional  parameter  space  and  3- 
dimensional  coefficient  space  will  be  used  to  illustrate  the  basic  concept  on  which  our 
algorithm  is  developed. 

In  Fig.5.2-l(a),  the  hypercube  3  is  a  square  with  edges  V2V3,  V3V4,  and 
V4Vj  parallel  to  either  8j -axis  or  82*axis  and  the  perturbed  parameter  domain  JS)  is  the 
straight  line  segment  connecting  V2  and  V4  In  Fig.5.2-l(b),  Xj,  X^  Xj,  and  X4  are  the 
images  of  the  vertices  Vj,  V2,  V3,  and  V4  of  3  respectively  in  the  3-dimensional 
coefficient  space.  The  polytope  IP (3),  i.e.,  the  convex  hull  of  the  four  image  points  Xr 
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X2,  X3,  and  X4,  is  a  pyramid  with  six  edges.  X2  and  X4  are  crucial  vertices  since  they  are 
in  <!(«©).  From  Theorem  5.2-3,  the  four  edges  XjX2,  X2X3,  X3X4,  and  X4Xj  are  the 
images  of  V  V2,  V2V3,  V3V4,  and  V4Vj  respectively  and  therefore  are  supplemental 
edges.  The  other  two  edges  of  the  polytope,  X}X3  and  X4X2,  are  fictitious  edges  which 
may  not  even  be  in  the  image  of  JS. 


y 


Fig. 5.2-1  Multilinear  mapping,  polytopes,  and  JS)j  partition  technique. 


The  image  of  £)  (J5>  resp.),  d( J9)  ( d(J9)  resp.),  can  be  constructed  as  follows.  Let 
Vs  and  Vb  be  the  center  points  of  the  line  segments  VjV4  and  V2V3  respectively.  It  is  easy 

to  see  that  the  line  segment  drawn  between  Va  and  Vfe,  i.e.,  V#Vb,  is  parallel  to  Vj  V2  and 
V4V3,  and  therefore  parallel  to  the  8j-axis.  Since  the  mapping  is  multilinear,  il(VtVb)  will 
be  the  straight  line  segment  X4Xb,  where  X#  =  and  v  =  ^(Vb)  are  the  center  points 
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of  the  line  segments  X}X4  and  X2X3  respectively.  Similarly,  $J(VxVy)  will  be  the  straight 
line  segment  X  X  ,  where  X  =  d(V  )  and  X  =  $(V  )  are  the  center  points  of  the  line 
segments  XjX2  and  X4X3  respectively.  It  is  easy  to  see  that  X#Xb  and  XxXy  are  both  in 
<J(JS)  and  the  intersection  point  of  XaXb  and  XxXy,  i.e.,  Xc  is  in  d(«0).  Repeating  the 

above  mapping,  we  can  see  that  d(JS))  and  £J ( J9 )  can  be  constructed  as  that  shown  in 

Fig.5.2-l(c).  In  Fig.5.2-l(c),  d(JS>)  is  a  saddle-shaped  surface  and  is  inside  the  pyramid 
XjX2X3X4,  i.e.,  1P(J§).  The  v  curve  in  d(JS>)  is  the  image  of  JS>. 

Note  that  d(J®)  c  d(J9)  c  &(£)).  There  is  an  easy  way  to  check  the  stability  of 
(P( J9)  and  therefore  a  sufficient  condition  for  the  stability  of  d(«©)  can  be  determined 
without  too  much  effort.  That  is,  d(J5))  is  stable  if  all  the  six  edges  of  TP(Js>):  XjX2, 
X2X3,  X3X4,  X4Xj,  XjX3  and  X4X2  are  stable  [52].  The  stability  of  the  line  segment 
XX  is  determined  by  using  the  fast  segment  stability  checking  algorithm  which  will  be 

summarized  in  Sec.5.4. 

In  the  proposed  algorithm  for  computing  the  stability  margin  kM,  we  need  to  check 
whether  a  subdomain  is  stable.  Consider  as  the  line  segment  V2V4  in  Fig.5.2-l(a) 
and  its  corresponding  polytope  as  shown  in  Fig.5.2-l(b).  J0.  is  stable  if  (P(^;)  is 

stable.  If  IPlJSj)  is  unstable  and  the  instability  is  caused  by  the  crucial  edges  or  vertices, 
then  J9.  is  unstable.  If  is  unstable  but  the  instability  is  not  caused  by  the  crucial 

edges  or  vertices,  then  the  information  is  not  enough  to  determine  the  stability  of  <©..  In 
this  case,  a  partition  technique  should  be  used  and  repeated  until  the  stability  of  JS);  is 

determined.  The  partition  technique  is  illustrated  briefly  as  follows. 

The  hypercube  J5).  (square  VjV2V3V4)  is  partitioned  equally  into  four 
hyperrectangles  (square  V4V#VcVy),  (square  VcVxV2Vb),  (square 
VaVjVxVc),  and  (square  VyVcVbV3),  by  two  hyperplanes  (line  VaVb  and  line  VxVy) 
which  are  orthogonal  to  the  axes  8j  and  62>  From  Fig.5.2-1,  it  is  easy  to  see  that  ^(JSK)  c 
^(JSj,)  u  *J(JSj2)  c  u  TPi&a)  where  the  polytopes  POis^)  and  flP(i)i2)  are  the 
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pyramids  X4X#XcXy  and  X^X^  respectively.  If  both  tPCJS^j)  and  are  stable 

then  40;  is  stable.  If  either  3P(J§.j)  or  75(40i2)  is  unstable  and  the  instability  is  caused  by 
the  crucial  edges  or  vertices  (note  that  X4,  X2,  and  Xc  are  a  crucial  vertices),  then  J9.  is 
unstable.  If  1P(^n)  or  1P(JS).2)  is  unstable  but  the  instability  is  not  caused  by  the  crucial 
vertices  or  edges,  then  the  JS^'s  corresponding  to  the  unstable  !P(jS).j)'s  should  be 
partitioned  further  and  the  partition  process  is  continued  until  the  stability  of  JS> .  is 

determined. 


Fig.5.2-2  Further  partitioning  in  the  parameter  space. 

From  Fig.5.2.2  we  can  see  that  the  perturbed  parameter  domain  J9.  (line  V4V2)  is 

inside  the  shaded  hyperrectangles  (squares)  and  as  the  number  of  the  shaded 
hyperrectangles  become  very  large  the  image  of  J0 .  will  be  close  to  the  union  of  the 

polytopes  corresponding  to  the  shaded  hyperrectangles.  In  our  algorithm,  the  number  of 
partitioning  usually  is  small  since  the  partition  process  is  terminated  whenever  all 

are  stable  or  some  crucial  vertex  or  edge  is  unstable.  The  partition  process  needs  to 
continue  only  if  some  flP(j§j.)'s  are  unstable  and  the  instability  is  not  caused  by  the  crucial 

vertices  or  edges.  Besides,  only  the  JS./s  related  to  the  unstable  F(J9.j)'s  need  to  be 
partitioned  further. 

In  the  next  section,  an  algorithm  based  on  the  above  theorems  for  computing  the 
real  MSM  kM  and  the  real  SSV  p  will  be  summarized. 
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5.3  Algorithm  For  The  Real  Structured  Singular  Value 

The  algorithm  is  developed  based  on  the  multilinear  mapping  between  the  parameter 
and  the  coefficient  spaces.  Some  recent  results  by  De  Gaston  and  Safonov  [48],  Sideris 
and  Pena  [50],  Bartlett,  Hollot,  and  Lin  [52],  and  Bouguerra,  Chang,  Yeh,  and  Banda  [57] 

are  used  in  the  algorithm.  The  iterative  algorithm  is  designed  to  find  the  largest  positive  real 
number  k,  i.e.,  the  real  MSM  kM,  such  that  the  system  in  Fig.5.1-1  remains  robustly  stable 

for  every  possible  perturbation  in  k^JS).  Then  the  real  SSV  jj.  can  be  obtained  from  |i  = 
1  /  kM- 

Recall  that  3  is  the  hypercube  with  2m  vertices  (±1, ....  ±1)  in  the  m-dimensional 
real  space  and  the  parameter  perturbation  domain  J®  is  an  r-dimensional  hyperplane  inside 
the  m-dimensional  hypercube  3.  In  Fig.5.1-1,  the  nominal  system  M(s)  is  fixed  and  the 
shape  of  the  perturbation  domain  J9  is  also  fixed.  In  the  computation  of  stability  margin, 
we  will  just  enlarge  or  shrink  the  perturbation  domain  J9  by  a  factor  k  until  k  is  the  largest 
positive  real  number  such  that  the  system  in  Fig.5.1-1  remains  robustly  stable  for  every 
possible  perturbation  in  kJ9.  Recall  that  "kJS)  is  stable"  means  "the  system  in  Fig.5.1-1 
remains  robustly  stable  for  every  possible  perturbation  in  k!9"  and  ”d(kJS>)  (resp.  1P(kj§)) 
is  stable"  means  "every  characteristic  polynomial  in  d(kJ9)  (resp.  !P(kJB))  is  stable." 

Since  d(kXD)  is  a  coefficient-space  image  of  kJ9,  kJS>  is  stable  if  and  only  if  d(kJ9) 
is  stable.  From  Theorems  2.6  and  2.7,  we  have  d(kJS>)  c  d(kj§) c  !P(kj0)  which  implies 
that  if  3P(kJ9)  is  stable  then  d(kJS)  is  stable  and  therefore  so  is  kJ®.  Hence,  any  k  such 
that  the  polytope  1P(kJ®)  remains  stable,  say  kL,  can  be  served  as  a  lower  bound  for  kM. 

The  stability  of  the  polytope  F(kJ9)  can  be  easily  checked  by  using  the  recent  results 
developed  by  Bartlett,  Hollot,  and  Lin  [52],  and  Bouguerra,  Chang,  Yeh,  and  Banda  [57]. 

Any  k  which  causes  instability  of  ki®  or  d(kJ®)  qualifies  as  an  upper  bound  for 
kM.  In  the  proposed  algorithm,  some  k  which  causes  instability  at  the  crucial  vertices  or 

edges  of  the  polytopes  corresponding  to  kj®  or  its  subdomains,  say  ky,  will  be  used  as  an 
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upper  bound  for  kM.  If  the  lower  and  the  upper  bounds  coincide  or  close  enough,  we  have 
the  stability  margin  kM.  Otherwise,  the  iterative  algorithm  needs  to  continue  until  the  gap  of 

the  lower  and  the  upper  bounds  is  small  enough. 

In  an  iteration  loop  of  the  algorithm,  for  a  given  k  we  construct  a  polytope  !P(kJ9). 
If  <P(ki9)  is  stable,  kL  is  updated  by  k.  In  the  case  of  unstable  3P( kjjj),  the  instability  may 

be  caused  by  the  crucial  or  the  noncrucial  edges.  If  it  is  caused  by  the  crucial  vertices  or 
edges  the  upper  bound  ky  shall  be  updated  by  k.  On  the  other  hand,  the  instability  is 

caused  by  the  noncrucial  edges,  no  information  can  be  used  to  update  kL  or  ky.  In  this 
situation,  the  perturbation  domain  partition  technique  [48,50]  is  employed. 


The  domain  partition  procedure  is  described  as  follows.  Refer  to  Fig.5.3-1,  the 
perturbation  domain  kJ9  is  the  line  segment  VaVd  and  kjj)  is  the  hypercube  (square)  which 
encloses  kJ9.  Assume  kL<©  (line  segment  VbVc)  is  stable  and  kjj)  is  unstable  and  the 

instability  is  caused  by  noncrucial  vertices  or  edges.  To  determine  the  stability  of  kJB,  the 
following  partition  technique  is  used.  The  domain  k  JS)  is  partitioned  into  three  parts:  kLJ9 

(line  segment  VbVc),  J$>j  (line  segment  V#Vb),  and  JS>2  (line  segment  VcVd).  Enclose  J0j 
and  J92  by  hyperrectangles  J9j  and  JS2  respectively.  Now,  there  are  three  possibilities:  (1) 
Both  VHJSj)  and  P(J92)  are  stable,  then  k JS)  is  stable  and  the  lower  bound  kL  shall  be 
updated  by  k.  (2)  Some  crucial  vertex  or  edge  of  V’fJDj)  or  P(JS2)  is  unstable,  then  k JS>  is 
unstable  and  the  upper  bound  ky  shall  be  updated  by  k.  (3)  Either  or  both  of  ^(JSj)  and 
iP(JS2)  is  unstable  and  the  instability  is  not  caused  by  crucial  vertices  or  edges,  then  no 
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information  can  be  used  to  update  kL  or  ky  and  therefore  the  JS/s  with  unstable  IP 

shall  be  partitioned  further  until  all  polytopes  which  contain  parts  of  d(kJS))  are  stable 
(update  the  lower  bound  kL  by  k)  or  some  crucial  vertex  or  edge  is  unstable  (update  the 

upper  bound  ky  by  k).  The  iterative  procedure  is  repeated  with  a  new  k  =  (ky+kL)/2  or  k  = 
2kL  (if  ky  is  not  available)  until  ky  -  kL  is  negligible  and  then  the  stability  margin  is  kM  = 

kL- 

The  algorithm  for  computing  the  real  MSM  kM  and  the  real  SSV  p  is  listed  as 
follows. 

ALGORITHM  5.3-1: 

1 .  Initialize  both  the  lower  and  the  upper  bounds  kL  =  ky  =  0.  Set  the  partition  indicator 
p  =  0.  Set  e  for  the  accuracy  of  k^  Set  initial  trial  value  at  k=l . 

2 .  Check  the  stability  of  the  corresponding  polytope  IP(kjs).  If  it  is  stable,  go  to  3. 
Otherwise,  go  to  5. 

3 .  Update  kL  by  k  and  check  ky.  If  ky  =  0,  update  k  by  2kL  and  go  to  2.  If  ky  *  0,  (it 
must  have  been  updated  and  must  ky  >  kL )  go  to  4. 

4 .  Set  k  =  (ky+kL)/2.  If  p  =  0,  then  go  to  2.  Otherwise  go  to  5. 

5 .  If  the  instability  of  the  polytope  !P(kJ9)  (or  the  polytopes  of  the  relevant  sub¬ 
hyperrectangles  as  resulted  from  the  partitioning  in  7)  is  due  to  crucial  vertices  or 
crucial  edges  (ki9  is  unstable),  go  to  8.  Otherwise,  go  to  6. 

6 .  Partition  the  hyperrectangles  with  unstable  polytopes  by  the  hyperplanes  orthogonal 
to  the  axes  of  coordinates.  Only  the  sub-hyperrectangles  which  enclose  kJ9  -  kLJS 

(we  call  these  the  relevant  sub-hyperrectangles)  need  to  be  considered. 

7 .  If  the  polytopes  of  the  relevant  sub-hyperrectangles  under  consideration  are  all  stable, 
set  p  =  1  and  go  to  3.  Otherwise,  go  to  5. 

8.  kJD  is  unstable.  Update  ky  by  k  and  check  ky  -  kL.  If  ky  -  kL  >  e,  go  to  4. 
Otherwise,  set  kM  =  kL  and  stop. 


68 


In  the  algorithm,  e,  the  tolerable  gap  between  the  lower  bound  kL  and  the  upper 
bound  kut  is  a  small  number  preset  to  determine  the  computation  accuracy.  The  algorithm 
is  employed  to  determine  the  stability  margin  kM  by  narrowing  down  the  gap  between  the 
lower  bound  kL  and  the  upper  bound  ky.  The  perturbation  domain  partitioning  loop  is 
described  in  step  6.  Only  the  hyperrectangles  enclosing  k«©  -  kLJS  and  corresponding  to 

unstable  polytopes  whose  crucial  vertices  and  edges  are  stable  need  to  be  partitioned  further 

and  the  domain  partitioning  loop  is  terminated  whenever  all  the  polytopes  under 
consideration  are  stable  or  any  crucial  vertex  or  edge  is  unstable.  The  set  kJ9  -  kLJS> 

consists  all  the  points  in  k J®  which  is  not  in  kLJ® .  Theoretically,  the  domain  partitioning 

loop  should  be  terminated  in  finite  steps  unless  the  unstable  region  is  of  measure  zero  (e.g. 
a  singular  point)  inside  k  JS) .  However,  in  practical  computation  a  preset  counter  can  be 
placed  in  the  loop  to  avoid  the  loop  from  repeating  too  many  times  for  a  given  k.  The 
partition  indicator  p  initially  is  set  as  zero  and  is  set  to  be  one  in  step  7  whenever  kJS)  is 
determined  to  be  stable  after  partitioning.  If  a  kLJ3  is  determined  to  be  stable  after 
partitioning,  then  the  polytope  <P(kJ9)  is  unstable  for  k  >  kL>  Therefore,  in  step  4  after 

setting  a  new  k  we  go  to  5  instead  of  2  if  p  *  0. 

Two  examples  are  used  to  illustrate  our  algorithm.  The  first  example  is  from  De 
Gaston  and  Safonov  [48]  in  which  J9  =  j® ,  i.e.,  the  elements  in  the  parameter  vector  8  are 
independent  and  the  other  example  has  dependency  among  the  entries  of  8.  In  the 

following,  the  vertices  of  the  hypercube  kj®  and  their  image  points  in  the  coefficient  space 
are  denoted  by  V/s  and  X.'s  respectively.  The  line  segment  between  X.  and  X.  is 
represented  by  X.X..  Note  that  both  V.'s  and  X;'s  are  functions  of  k. 


Fig.  5.3-2  A  closed-loop  system  considered  in  Example  5.3-1. 
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Fig.  5.3-3  Perturbation  domain  of  Example  5.3-1 


Example  5.3-1  : 

De  Gaston  and  Safonov  [48]  developed  an  algorithm  for  computing  the  stability 
margin  based  on  the  multilinear  mapping  between  the  parameter  space  and  the  complex 
plane  [61].  They  used  an  example  to  illustrate  their  algorithm.  We  will  use  exactly  the  same 
example  to  illustrate  our  algorithm  and  then  the  solutions  of  both  approaches  can  be 
compared. 


The  system  to  be  considered  is  shown  in  Fig.5.3-2  where  the  plant  P(s)  and  the 
controller  C(s)  are  described  by 


P(s)  = 


i« 


s  ( s  +  82a)  (s  +  83i) 


and 


C(s)  = 


s  +  z 
s  +  p 


i_ 

l 


respectively.  The  parameters  in  the  above  expressions  arc  given  by 


Zj=  2  rad/sec 

Pj=  10  rad/sec 

5u=8,o(l+8i) 

On 

o 

II 

00 

8 

1 8, 1  £  0.1 

S2.=  52o+  52 

8^  =  4  rad/sec 

1 82 1  0.2 

rad/sec 

S3.=  83o+  53 

8^=  6  rad/sec 

1 83 1  £  0.3 

rad/sec 

The  perturbation  part  of  the  perturbed  closed-loop  system  can  be  pulled  out  and  represented 
by  a  diagonal  matrix  ar.d  therefore  the  perturbed  closed-loop  system  in  Fig.  5.3-2  can  be 


restructured  as  that  in  Fig. 5. 1-1,  where  A  =  diag  {  8,,82,  83  }  and  the  nominal  closed- 
loop  system  M(s)  has  a  realization  (A,  B,  C)  as  follows: 
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The  vertices  of  the  perturbation  domain  k  JS>  in  the  3-dimensional  parameter  space 
are  numbered  as  shown  in  Fig. 5. 3-3.  Recall  that  the  vertices  of  kJS)  are  denoted  by  V,'s. 

The  vertices  of  k  JS)  with  k  =  1  are 


V,=  (-.l,-.2,  -.3), 
V,=  (.1,  -.2,  .3), 


V2=  (-.1,  .2,  -.3), 
V6=  (.1,  .2,  .3), 


V3=  (-.1,  .2,  .3), 
V7=  (.1,  .2,  -.3), 


V4=  (-.1,  -.2,  .3), 
Vg=  (.1,  -.2,  -.3). 


Since  M(s)  is  a  rational  matrix  with  real  coefficients,  the  coefficients  of  the  characteristic 
polynomial  are  multilinear  functions  of  the  parameters  8,,S2,  and6y  The  characteristic 

polynomial 

aQs4  +  a,s3  +  a2s2  +  a3s  +  a4 

can  be  obtained  from  det  (si  -  (A  -  BAC)]  and  then 
<x0  =  1 

a,  =  20  +  82  +  83 

a2  =  124  +  1682  +  1483  +  8283 

a3  =  1040  +  8008,  +  6082  +  4083  -t-  10S283 

a4  =  1600  +  16008, 

As  we  expect,  the  coefficient  vector  a  =  (aQ  a,  a4]T  is  a  multilinear  function  of  the 
parameter  vector  8  =  (8,  82  83]T. 


The  objective  is  to  compute  the  stability  margin  kM  which  is  the  largest  k  such  that 

kJS)  is  stable,  i.e.,  the  closed-loop  system  remains  robustly  stable  in  kJ9.  For  each  k,  we 
map  the  eight  vertices  of  kJB  into  the  coefficient  space  as  X.,  i=l,2,...,8  from  which  the 

polytope  IP(kJS))  can  be  constructed.  It  is  easy  to  see  that  the  polytope  IP(kJS))  has  12 
crucial  edges  and  16  fictitious  edges.  If  these  28  edges  are  all  stable,  i.e.,  3P(kJ9)  is  stable, 
then  kJS  is  stable.  If  any  of  the  crucial  edges  is  unstable,  then  so  is  kJ9.  Note  that  unstable 
lP(kJ0)  does  not  imply  unstable  kJS,  since  the  instability  may  be  only  caused  by  some 
fictitious  edges  which  are  not  in  tJ(kJ0). 

Let's  go  through  the  algorithm  listed  in  the  previous  section.  First  of  all,  we 
initialize  both  the  lower  and  the  upper  bounds  kL  =  ky  =  0  and  set  the  accuracy  e  =  10‘6. 

For  k=l,  1P(kJS>)  is  stable  and  we  have  kL  =  1  and  update  k  to  2.  For  k=2,  again  1P(kJ0)  is 
stable  and  therefore  we  update  kL  to  2  and  k  to  4  respectively.  For  k=4,  P(kJ9)  is  unstable 
and  the  instability  occurs  at  Xg  which  is  crucial  and  ky  is  updated  to  4. 

Now,  we  are  in  the  loop  of  bisection  to  narrow  down  the  gap  between  kL  and  ky  to 
be  less  than  e.  After  21  iterations  we  have  kL  =  3.417395  and  ky  =  3.417396  such  that 
HP(kL«C))  is  stable  while  IPfkyJS))  is  unstable  and  the  instability  occurs  at  Xg.  Therefore,  we 
have  the  real  MSM  kM  =  3.417395. 

From  the  above  computation,  we  find  that  the  no.  8  vertex  of  the  perturbation 
domain  is  the  critical  point.  Thus,  we  can  check  our  solution  at  Vg  of  the  parameter  space. 

When  k  =  kM  =  3.417395,  the  closed-loop  system  at  Vg  is  stable  since  its  characteristic 

values  are 

-16.3521921088081 

-.00000016598452  +  8.22820065594421 
-.00000016598452  -  8.22820065594421 
-1.93911005922287 

If  we  increase  k  a  little  bit  to  k  =  ky  =  3.417396,  the  closed-loop  system  at  Vg  becomes 
unstable  since  its  characteristic  values  are 
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-16.352192213492 

.00000014646442  +  8.2282008933441 
.00000014646442  -  8.2282008933441 
-1.9391100794366 


In  [48],  De  Gaston  and  Safonov  obtain  a  margin  k  =  3.44  which  is  close  to  our  result. 
However,  when  k  =  3.44,  the  closed-loop  system  at  Vg  is  unstable  since  its  characteristic 


values  are 


-16.354555659044 

.00706068635988  +  8.23356355623191 
.00706068635988  -  8.23356355623191 
-1.9395657136757 


That  means  3.44  could  not  be  a  correct  margin.  One  of  the  major  reasons  that  De  Gaston 
and  Safonov's  algorithm  [48]  is  complicated  is  that  it  requires  frequency  ©  sweep.  For  the 
same  example,  they  found  that  at  ©  =  8.22  rad/sec,  det[I+kM(j©)A]  approximately  equals 
to  zero  at  Vg  when  k  =3.44.  To  have  an  accuracy  of  10"6  for  the  stability  margin,  they  need 

a  more  accurate  ©,  say  ©  =  8.2282  rad/sec,  which  needs  extensive  frequency  ©  search. 
Our  approach  does  not  require  frequency  ©  search. 

The  real  SSV  |l  for  the  system  is|X=l/kM=l  /  3.417395  =  .29262055. 


Example  5.3-2  : 

Assume  that  the  nominal  system  M(s)  in  Fig.5.1-1  has  a  state-space  representation 
(A,B,C)  as 
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and  the  perturbation  matrix  A  in  Fig.5.1-1  is  given  by 

A  =  diag  {  5j,  52  ,  63  } 

where 
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and 


8 


1 


-liSSjSl,  -1  £  82£  1 

The  set  of  the  characteristic  polynomials  of  the  perturbed  system  can  be  described  by 

aQs4  +  (XjS3  +  a2s2  +  a3s  +  a4 


where 
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Fig.5.3-4  Perturbation  domain  of  Example  5.3-2. 


The  parameter  perturbation  domain  is  shown  in  Fig.5.3-4.  The  shaded  area  V,V0V,V,  js 
the  perturbation  domain  J9  in  which  8j  =  8 1 ,  82  ■  83  *  82  and  the  hypercube 
V!V2V3V4V5V6V7V8  is  ^ 1116  following  we  will  use  the  proposed  iterative  algorithm 
to  compute  kM.  Initially,  kL  =  ky  =  0  and  e  is  set  as  104.  P(kJS)  is  stable  for  k=l  and 
therefore  kL  is  updated  to  1.  For  k=5,  the  polytope  P(kJS))  is  unstable  and  the  instability  is 
caused  by  the  crucial  vertices  X3,  Xfi  and  by  the  noncrucial  vertex  X?. 
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Therefore  ky  is  updated  to  5.  As  a  next  step  k  is  chosen  as  (l+5)/2  =  3.  FfkjS)  is 
stable  for  k=3  and  therefore  kL  is  updated  to  3.  Following  the  iteration  steps,  k=(3+5)/2=4, 

we  find  that  the  polytope  IP(kJS))  is  unstable  for  the  value  of  k  =  4  and  the  instability  is 
caused  by  the  crucial  vertex  Xg.  Therefore  ky  is  updated  to  4.  Iterating  k  between  3  and  4 

as  it  is  suggested  in  steps  2-8  we  find  that  the  polytope  <P(kJS))  is  unstable  for  the  value  of 
k  =  3.6297  and  the  instability  is  caused  by  the  crucial  vertex  Xft.  Hence,  ky  =  3.6297. 

IP(kJS))  is  stable  at  for  k=3.6296  and  we  have  kL=  3.6296.  Now,  Iky-  kj  <,  e  and 
therefore  we  have  the  real  MSM  kM  =  3.6296.  The  real  SSV  isji=l/kM  =  l/  3.6296  = 
0.2755  . 

5.4  Fast  Algorithm  to  Check  the  Stability  of  Segments 

In  the  previous  subsection,  an  iterative  algorithm  for  computing  the  real  SSV  based 
on  poly  topic-polynomial  approach  was  presented.  In  this  algorithm,  the  major  computation 
is  checking  the  stability  of  edges  of  polytopes. 

There  are  several  different  computational  methods  available  to  verify  the  stability  of 
an  exposed  edge  of  a  polytope.  However,  in  practical  problems,  these  methods  are  found 
to  impose  a  heavy  computational  burden  due  to  the  fact  of  the  "combinatoric  explosion"  of 
the  number  of  edges  of  the  polytope  whose  stability  needs  to  be  checked.  To  illustrate  this 
phenomenon,  we  consider  a  system  with  m  perturbation  parameters.  The  corresponding 
polytope  in  the  coefficient  space  has  2m  vertices.  There  are  2m  l(2m-l)  edges  joining  these 
vertices;  note  that  the  number  of  edges  which  depend  on  m  can  be  very  large  compared  to 
the  number  of  vertices.  This  is  a  fact  that  will  be  exploited  later  on  in  our  algorithm  in  order 
to  reduce  the  computations  to  a  minimum. 

These  existing  methods  to  verify  the  stability  of  these  edges  usually  require 
laborious  computations  that  need  to  be  performed  independently  for  each  edge.  For 
instance,  there  is  the  roots  locus  technique  which  consists  of  i  presenting  the  edge  of  a 
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polytope  as  a  convex  combination  of  two  polynomials  a(s)  and  (5(s)  of  degree  n.  The 
characteristic  polynomial  of  this  convex  combination  can  be  written  as  1  +  k(3(s)/a(s), 
where  k  takes  on  the  values  from  zero  to  infinity.  The  necessity  to  sweep  the  infinite 
interval  for  k  makes  this  method  unattractable.  Another  widely  used  test  for  the  stability  of 
the  convex  combination  of  two  stable  polynomials  was  developed  by  Bialas  [54].  He 
showed  that  the  convex  combination  of  two  stable  polynomials  of  degree  n  is  stable  if  and 
only  if  all  the  real  eigenvalues  of  -H(a)H'1([3)  are  negative,  where  H(a)  is  the  Hurwitz 
matrix  associated  with  the  polynomial  a(s).  Although  Bialas'  test  does  not  require  the 
sweep  of  k,  the  computations  involved,  i.e.  computations  of  eigenvalues  of  an  nxn  matrix, 
must  be  carried  out  for  each  edge  of  the  polytope  independently.  Hence,  in  the  case  of 
combinatoric  explosion  alluded  to  earlier,  this  technique  can  be  found  to  be  computer  time 
consuming.  The  iterative  algorithm  proposed  by  [60]  for  computing  the  multivariable 
stability  margin  needs  to  check  the  stability  of  polytopes  for  each  iteration.  Therefore,  a 
more  efficient  computational  tool  for  checking  the  stability  of  the  edges  of  a  polytope  is 
necessary. 

In  this  section,  we  propose  a  fast  algorithm  for  checking  the  stability  of  the  edges  of 
a  poly  tope  where  most  of  the  computations  involved  depend  on  the  number  of  vertices 
rather  than  on  the  number  of  edges.  This  algorithm  was  developed  based  on  "the  segment 
lemma"  derived  by  Chapellat  and  Bhattacharyya  U>8].  Although  the  segment  lemma  is  a 
great  result,  no  explicit  algorithm  was  given  in  [58].  We  will  reveal  some  magnificent 
properties  of  the  lemma  and  show  how  these  will  lead  to  a  fast  algorithm.  In  the  proposed 
algorithm,  the  major  computation  involved  is  the  solution  of  the  positive  real  roots  of  two 
polynomials  with  degree  less  than  or  equal  to  n/2  for  each  vertex.  The  computation  required 
by  the  algorithm  is  mainly  vertex-dependent,  and  the  burden  of  the  combinatoric  explosion 
of  the  number  of  edges  is  greatly  reduced. 
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The  major  computations  in  the  algorithm  are  vertex  dependent,  and  the  edge 
dependent  computations  involved  are  trivial.  Therefore  the  combinatoric  explosion  is  not 
fatal  if  the  proposed  algorithm  is  used. 

As  mentioned  earlier,  the  problem  of  checking  robust  stability  of  a  system  can  be 
reduced  to  that  of  checking  the  stability  of  the  edges  of  polytopes.  However,  an  edge  of  a 
polytope  is  nothing  but  a  line  segment  in  the  coefficient  space  which  can  be  represented  as 
the  convex  combination  of  two  vertices  corresponding  to  two  polynomials.  If  the  two  end 
points  of  the  line  segment  are  designated  by  a  and  |3,  hence  corresponding  to  two 
polynomials  a(s)  and  (3(s)  of  degree  n,  then  every  point  on  the  line  segment  corresponds  to 
a  polynomial  is  given  by 

t  a(s)  +  (l-t)P(s)  (5.4-1) 

for  some  t  e  [0,1].  In  the  following,  a(s)  and  j3(s)  are  assumed  stable.  Let  the  polynomials 
a(s)  and  p(s)  of  degree  n  be  written  as 

a(s)  =  a(s)  +  s  b(s)  (5.4-2a) 

P(S)  =  c(s)  +  s  d(s)  (5.4-2b) 

where  the  polynomials  a(s),  b(s),  c(s),  d(s)  are  of  the  form 

p(s)  =  pQ  +  pj  s2  +  p2  s4  +  p3  s6  +  -  (5.4-3) 

Instability  of  the  line  segment  occurs  when  the  roots  of  (5.4-1)  cross  the  imaginary  axis. 
Therefore,  letting  s  =  jco  and  substituting  x  for  ci)2,  (5.4-3)  becomes 

£(x)  =  p0  -  Pj  x  +  p2  x2  -p3  x3  +  -  (5.4-4) 

where  the  degree  of  p(x)  is  now  only  half  of  that  of  p(s).  Setting  (5.4-1)  equal  to  zero, 
substituting  (5.4-2)  and  (5.4-4)  and  equating  real  and  imaginary  parts  on  both  sides  of  the 
resulting  equation,  yields 
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and  therefore 


t  6(x) 

1-t  S(x) 

(5.4-5a) 

t  _  to 

i-‘  '  6(x) 

(5.4-5b) 

fi(x)  3(x)  -  6(x)  c(x)  =  0 

(5.4-6) 

where,  for  n  even,  the  polynomials  a(x),  6(x),  6(x),  and  3(x)  are  expressed  as: 


a(x)  = 

S(x)  = 


8(x)  = 


^  (-l)kc.xk 
k=0  k 


s 

k=0 


(-  i)kvk 


(5.4-7) 


(5.4-8) 


(5.4-9) 


(5.4-10) 


when  n  is  odd  the  summation  in  the  above  four  equations  is  carried  out  up  to  (n-l)/2. 
Equation  (5.4-6)  then  becomes 


L  i(-l)k(»iVi-bl,ci)  x“  =  0 
k=0  i=0 


note  that  equation  (5.4-1 1)  is  of  order  n-1. 


(5.4-11) 


Now,  according  to  the  segment  lemma  [58],  the  line  segment  (edge)  connecting  the 
two  stable  points  a  and  (3  in  the  coefficient  space  is  unstable  if  and  only  if  there  exists  a 
positive  real  x  which  satisfies  (5.4-6)  and  (5.4-5)  simultaneously.  One  can  solve  the 
positive  real  roots  of  the  (n-l)-th  order  polynomial  equation  (5.4-6)  first  and  then  plug 
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these  roots  into  either  (5.4-5a)  or  (5.4-5b).  If  one  of  these  roots  make  the  left  handed  side 
of  (5.4-5),  i.e.,  t/(l-t),  positive,  then  the  edge  is  unstable.  Otherwise,  the  edge  is  stable. 
This  approach  is  quite  straightforward  but  it  is  still  not  efficient. 

Recall  that  by  using  Bialas'  method  [54],  one  has  to  construct  a  Hurwitz  matrix  H. 
of  size  n  and  invert  the  matrix  for  each  vertex,  then  compute  the  product  H.H.'1  and  find 

the  eigenvalues  of  the  nxn  matrix  for  each  edge.  It  is  easy  to  see  that  the  approach  we 
mentioned  in  the  previous  paragraph  is  slightly  better  than  Bialas'  since  solving  positive 
real  roots  of  an  (n-l)-th  order  polynomial  equation  is  easier  than  computing  the  eigenvalues 
of  an  nxn  matrix.  However,  both  approaches  require  nontrivial  computations  for  each  edge 
of  the  polytope  and  hence  they  will  become  impractical  when  the  combinatoric  explosion  of 
the  number  of  edges  occurs. 

In  the  proposed  algorithm,  first  of  all,  we  solve  the  positive  real  roots  for  the  four 
polynomials  a(x),  6(x),  c(x),  and  3(x),  Then  from  these  positive  roots,  the  intervals  such 
that  both  a(x)c(x)  and  cl(x)6(x)  are  negative  can  be  easily  found.  If  such  intervals  do  exist 
one  needs  to  check  for  the  existence  of  roots  of  equation  (5.4-6)  inside  these  intervals.  If 
(5.4-6)  admits  roots  inside  these  intervals  then  the  edge  is  unstable  and  stable  otherwise. 
Note  that  when  such  intervals  are  empty  one  needs  not  continue  and  can  conclude  that  the 
edge  is  stable.  Recall  that  when  n  is  even  a(x)  and  c(x)  are  of  order  n/2  and  6(x)  and  cl(x) 
are  of  order  n/2-1.  When  n  is  odd,  a(x),  6(x),  6(x),  and  S(x)  are  all  of  order  (n-l)/2.  One  of 
the  significant  features  of  this  method  is  revealed  in  the  fact  that  the  polynomials  a(x),  6(x), 
c(x),  and  cl(x)  are  of  at  most  degree  n/2,  hence  it  is  computationally  easier  to  deal  with.  In 
fact,  if  n=4,  a(x)  and  c(x)  are  second-order  polynomials  and  £>(x)  and  S(x)  are  first-order 
polynomials.  If  n=9,  a(x),  c(x),  6(x),  and  3(x)  are  fourth-order  polynomials.  From  [63], 
we  can  have  closed-form  solutions  for  the  polynomial  equations  with  order  less  than  or 
equal  to  four.  Also,  the  problem  of  determining  whether  or  not  (5.4-6)  has  positive  real 
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roots  inside  an  interval  is  a  trivial  one  since  it  can  be  easily  checked  by  using  the  Sturm 
Theorem  [see  App.  B]  without  solving  the  equation. 

The  most  significant  feature  of  the  proposed  algorithm  is  that  the  computation 
involved  is  mainly  determined  by  the  number  of  vertices  instead  of  that  of  edges.  To  see 
why  the  algorithm  is  more  efficient  than  Bialas',  let  us  assume  the  number  of  the 
perturbation  parameters  m=10  and  the  order  of  the  characteristic  polynomials  n=9.  Now  in 
the  coefficient  space,  the  polytope  has  210=1024  vertices  and  more  than  a  half  million 
(523,776)  edges.  By  using  Bialas'  method  to  check  the  stability  of  the  polytope,  we  need 
to  construct  a  9x9  Hurwitz  matrices  R  and  its  inverse  H.'1  for  every  vertex  and  then  do  the 
matrix  multiplication  HR1  and  find  the  eigenvalues  of  the  9x9  matrix  HR'1  for  each 
edge.  By  using  the  new  algorithm,  for  each  vertex  we  only  need  to  solve  two  4-th  order 
polynomials  for  positive  real  roots  and  then  for  each  edge  find  the  intervals  such  that  both 
a(x)c(x)  and  S(x)6(x)  are  negative  which  can  be  easily  found.  If  such  intervals  are  empty, 
the  edge  is  stable.  On  the  other  hand,  if  such  intervals  do  exist  we  need  to  use  the  Sturm 
theorem  to  check  if  the  8-th  order  equation  (5.4-11)  has  any  roots  inside  these  intervals. 
Now  we  can  see  that  the  major  computation  effort  in  the  new  algorithm  is  finding  the 
positive  real  roots  for  2,048  4-th  order  polynomial  equations  while  Bialas'  method  requires 
finding  the  eigenvalues  for  more  than  a  half  million  9-th  order  matrices. 

Three  examples  are  used  to  illustrate  the  segment  stability  checking  algorithm. 
Example  5.4*1: 

Consider  the  two  stable  polynomials  given  as: 

a(s)  =  s4  +  5s3  +  3s2  +  2s+l 

J3(s)  =  s4  +  s3+5s2+  s  +  3 
The  polynomials  in  equations  (5.4-7)  -  (5.4-10)  become 

a(x)  =  1  -  3  x  +  x2 
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6(x)  =  2  -  5  x 
8(x)  =  3  -  5  x  +  x2 
S(x)  =  1  -  x 

Upon  substitution  into  equation  (5.4-1 1),  one  obtains 

4  x3  -  23  x2  -i-  21  x  -  5  =  0 

Instability  of  the  convex  combination  occurs  whenever  there  are  positive  values  of  x  that 
satisfy  both  equations  (5.4-5a)  and  (5.4-5b)  as  well  equation  (5.4-6).  For  this  example, 
one  can  see  that  equations  (5.4-5)  are  satisfied,  i.e.,  both  a(x)c(x)  and  cl(x)fe(x)  are 
negative  over  the  interval  [.4,  .697].  To  see  whether  (5.4-6)  admits  any  roots  inside  this 
interval,  the  following  Sturm  sequence  is  formed: 

f0  (x)  =  4  x3  -  23  x2  +  21  x  -  5 
fj  (x)  =  12  x2  -  46  x  +  21 
f2  (x)  =  554  x  -  303 
f3  (x)  =  87354 

Therefore  the  number  of  real  roots  of  (5.4-6)  inside  [.4,  .697]  is  the  difference  in  sign 
variation  of  the  sequence  when  evaluated  at  both  end  points,  i.e.  in  this  example  there  are 
two  roots  inside  the  interval  and  so  the  edge  is  unstable. 

Example  5.4-2: 

As  a  second  example  we  now  consider  the  following  two  polynomials: 

a(s)  =  s4  +  5  s3  +  10  s2  +  5  s  +  1 
(3(s)  =  s4  +  2  s3  +  15  s2  +  s  +  3 

The  polynomials  in  equations  (5.4-7)  -  (5.4-10)  become 

&(x)  =  1  - 10  x  +  x2 
6(x)  =  5  -  5  x 
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6(x)  =  3  -  15  x  +  x2 


<}(x)  =  1  -  2  x 

Upon  substitution  into  equations  (5.4-5),  one  can  see  that  a(x)c(x)  is  negative  over  [.101, 
.203]  and  [9.9,  14.79]  whereas  9(x)S(x)  is  negative  over  [.5, 1].  Hence  there  are  no  values 
of  x  that  satisfy  both  equations  (5.4-5a)  and  (5.4-5b)  at  the  same  time  and  therefore  the 
edge  is  stable.  Note  that  checking  for  the  existence  of  roots  of  equation  (5.4-6)  is  not 
required  for  this  example. 

Example  5.4-3: 

In  this  third  example,  we  apply  our  algorithm  to  the  perturbed  plant  in  Example  5.3-1.  In 
Example  5.3-1,  we  found  that  for  k=3.417396  the  characteristic  polynomial  corresponding 
to  Vg  is  unstable.  For  k=3.417395  the  polytope  corresponding  to  kJS)  is  stable.  This 

polytope  has  8  vertices  and  28  edges.  Recall  that  if  we  use  Bialas'  method  to  check  the 
stability  of  the  polytope,  we  need  to  construct  a  Hurwitz  matrix  H.  and  its  inverse  H.1  for 
each  vertex,  do  the  multiplication  RH:1  and  find  the  eigenvalues  of  a  4x4  matrix  H.R  1 

for  each  edge.  In  the  following,  we  will  use  our  fast  algorithm  to  verify  that  this  polytope  is 
stable.  We  first  solve  a  first-order  and  a  second-order  equations  associated  with  each 
vertex.  Then,  for  each  pair  of  vertices  (each  edge)  the  intervals  of  x  such  that  (5.4-5)  are 
satisfied  can  be  easily  found.  It  turns  out  that  in  this  example  such  intervals  are  all  empty 
and  so  all  the  edges  of  the  polytope  are  stable.  The  stability  of  kJS  is  then  inferred  by  the 
edge  theorem.  Note  that  since  no  intervals  satisfying  equations  (5.4-5)  were  found,  one 
needs  not  bother  with  checking  the  conditions  required  by  equation  (5.4-6).  With  this,  we 
conclude  that  the  new  approach  can  be  used  effectively  to  greatly  reduce  the  amount  of 
computations  required  especially  when  one  is  faced  with  a  large  number  of  perturbation 
parameters  in  the  plant 
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SECTION  6 


CONCLUSION  AND  FURTHER  RESEARCH 
6.1  Concluding  Remarks 

As  demonstrated  in  Section  2,  most  robust  control  problems  can  be  formulated  as 
the  standard  H“  optimization  problem.  The  standard  H“  optimization  problem  consists  of 
two  subproblems.  One  is  stabilization  and  the  other  is  H°°  optimization.  The  observer-based 
controller  parametrization  is  employed  to  characterize  the  set  of  all  stabilizing  controllers  in 
terms  of  a  proper  stable  parameter  matrix.  Then  among  this  set  of  stabilizing  controllers, 
some  will  be  chosen  to  solve  the  H°°  optimization. 

The  poles  of  the  closed-loop  system  with  the  observer-based  controller 
parametrization  can  be  classified  into  three  groups,  and  each  group  of  poles  can  be 
independently  determined.  These  three  groups  of  poles  are  the  regulator  poles,  the  observer 
poles,  and  the  poles  of  the  added  parameter  matrix.  The  regulator  gain,  the  observer  gain, 
and  the  added  parameter  matrix  are  free  parameters  to  be  chosen  such  that  the  closed-loop 
transfer  function  matrix  has  some  optimal  performance  with  the  constraint  that  the  regulator 
poles,  the  observer  poles,  and  the  poles  of  the  parameter  matrix  are  stable.  Furthermore, 
these  parameters  can  be  chosen  such  that  the  controller  is  uncontrollable  and/or 
unobservable  and  then  the  uncontrollable  and/or  unobservable  controller  poles  can  be 
removed  and  the  order  of  the  controller  can  be  reduced.  The  set  of  these  removable 
controller  poles  is  a  subset  of  the  regulator  and  the  observer  poles.  The  poles  of  the  closed- 
loop  system  with  the  minimal  order  controller  will  include  all  the  poles  of  the  parameter 
matrix  and  some  of  the  regulator  and  the  observer  poles  which  are  not  the  removable 
controller  poles. 

There  is  very  few,  if  there  is  any,  one-block  problems.  Most  of  practical  control 
problems  are  either  two-block  or  four-block  problems.  To  reduce  the  order  of  the  optimal 
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controller,  we  require  low  order  realizations  of  Rn(s),  R12(s),  R21(s),  and  R22(s)  instead 
of  a  minimal  realization  of  R(s).  In  Sec.4.2,  a  numerically  reliable  algorithm  for 
constructing  the  minimal  realizations  of  RLi(s),  Rl2(s)>  Rri(s),  «nd  Rr2(s)  has  been 
presented.  From  these  realizations,  the  rational  matrices  Rij(s)  =  Ru(s)RRj(s),  ij=l,2,  can 
be  easily  obtained.  There  is  no  possible  mathematically  identical  pole-zero  cancellation 
between  RLi(s)  and  Rrj(s),  therefore,  the  order  of  the  realization  of  Rjj(s)  constructed  by 
our  algorithm  is  minimal  in  most  practical  cases. 

The  two-block  H°°  optimization  problem  can  be  solved  by  using  the  proposed  fast 
y-iteration  algorithm  in  Sec.4.3  to  compute  the  optimal  H°°  norm.  The  fast  y-iteration 
algorithm  converges  extremely  fast.  The  optimal  H°°  norm  can  be  obtained  with  accuracy 
up  to  double  precision  within  four  iterations.  Once  the  optimal  H“  norm  is  available,  the 
two-block  H“  optimization  problem  can  be  converted  into  a  one-block  H°°  optimization 
problem  and  then  an  optimal  controller  can  be  easily  constructed  by  using  Glover's 
method. 

For  the  four-block  H“  optimization  problem,  the  easiest  solution  is  probably  the 
two-Riccati-equation  approach  proposed  by  Doyle,  Glover,  Khargoneckar,  and  Fiancis. 
Strictly  speaking,  the  two-Riccati-equation  approach  only  gives  suboptimal  controllers  such 
that  the  H~  norm  of  the  closed-loop  system  is  less  than  a  number  y  which  is  larger  than  the 
optimal  norm.  In  Sec.4.4,  we  attempted  to  use  the  two-Riccati-equation  approach  to  design 
a  nearly  H“  optimal  controller.  A  simple  iterative  scheme  can  be  used  to  reduce  y  to  a 
number  which  is  very  close  to  the  optimum.  However,  as  y  is  close  to  the  optimum,  the 
elements  of  the  state-space  realization  of  the  controller  will  approach  to  infinity.  The 
remedy  to  this  difficulty  is  to  do  partial  fraction  expansion  for  the  controller  and  to 
approximate  the  wide-band  low-pass  terms  by  direct  feedthrough  terms.  Then  the  elements 
of  the  realization  of  the  controller  will  be  all  finite,  and  the  order  of  the  nearly  H°°  optimal 
controller  is  at  least  one  less  than  that  of  the  plant. 
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The  real  structured  singular  value  computation  algorithm  presented  in  Sec.5.3  was 
developed  based  on  the  multilinear  mapping,  the  edge  theorem,  the  fast  segment  stability 
checking  algorithm  in  Sec.5.4,  and  the  perturbation  domain  partition  technique.  The 
computational  burden  due  to  the  combinatoric  explosion  of  the  number  of  edges  has  been 
reduced  a  lot  by  using  the  fast  segment  stability  checking  algorithm.  The  convergence 
usually  is  not  a  problem  unless  a  singular  unstable  region  is  inside  or  very  close  to  the 
perturbation  domain  which  hardly  happens  in  the  real  world.  If  that  happens,  we  just 
simply  shrink  the  perturbation  domain  and  try  to  find  the  largest  lower  bound  kL  at  which 

the  computer  can  handle  in  finite  steps.  Although  the  algorithm  requires  no  frequency 
search,  the  frequency  at  which  the  system  becomes  unstable  can  be  computed  from  the  real 
structured  singular  value  and  the  critical  point  at  which  the  stable  perturbation  domain 
touches  the  unstable  region.  This  critical  point  can  be  evaluated  by  using  the  fast  segment 
stability  checking  algorithm. 

6.2  Further  Research  Suggestions 

The  two-Riccati-equation  approach  is  a  great  break-through  in  the  solution  of  H°° 
optimization  problem.  However,  it  only  gives  suboptimal  solutions.  If  we  try  to  have  a 
suboptimal  solution  which  is  very  close  to  the  optimum,  an  inevitable  numerical  difficulty 
will  occur.  This  numerical  difficulty  is  caused  by  an  unnecessary  restriction  that  the 
controller  be  strictly  proper.  The  first  further  research  suggestion  is  to  develop  two-Riccati- 
equation  like  state-space  formulae  for  optimal  (instead  of  suboptimal)  H“  controllers  which 
are  allowed  to  have  direct  feedthrough  terms. 

The  second  further  research  suggestion  is  related  to  the  M-A  structure  shown  in 
Fig.5.1-1.  M-A  structure  is  important  since  all  the  SSV  analysis  and  design  techniques  are 
developed  based  on  this  structure.  The  construction  of  an  M-A  structure  for  a  perturbed 
system  is  not  difficult  at  all.  However,  it  is  by  no  means  an  easy  task  to  find  a  minimal  M- 
A  structure.  Here  "minimal"  means  that  the  dimension  of  A  (or  M)  is  minimal.  We  can 
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easily  construct  an  M-A  structure  for  a  given  perturbed  system,  but  the  dimension  of  the 
structure  may  be  unnecessarily  large.  Although  theoretically  the  structured  singular  values 
of  all  M-A  structures  of  the  same  perturbed  system  should  be  identical,  a  minimal  M-A 
structure  is  much  easier  to  handle  than  those  with  high  dimension  in  computation.  A 
minimal  M-A  structure  is  essential  in  the  computation  of  the  structured  singular  value. 

The  third  further  research  suggestion  is  the  development  of  a  controller  design 
updating  scheme.  The  H°°  optimization  technique,  fast  computational  algorithm  for  the  real 
SSV,  and  the  controller  updating  procedure  are  employed  in  the  proposed  design  approach 
for  robust  optimal  controllers.  The  H°°  optimization  technique  allows  us  to  consider  a  larger 
and  more  realistic  set  of  disturbances,  noises,  and  commands  than  those  considered  by 
LQG  and  Wiener-Hopf  methods.  The  robust- stability  constraint  for  the  unstructured  plant 
uncertainty  and  the  control-input  constraint  are  also  incorporated  in  the  H”  cost  function. 
The  controller  which  minimizes  the  H°°  cost  function  is  a  tentative  robust  optimal  controller. 
For  the  closed-loop  system  with  the  tentative  robust  optimal  controller,  we  will  check  all 
performance  measures  like  the  real  SSV,  the  H°°  norm  of  the  sensitivity  function,  the  H°° 
norm  of  the  complementary  sensitivity  function,  and  the  amplitude  of  the  control  input,  etc. 
These  performance  measures  will  be  used  to  update  the  H~  cost  function  (and  therefore  the 
controller)  and  the  updating  procedure  is  repeated  until  a  satisfactory  trade  off  is  reached. 
In  the  proposed  design  approach,  the  uncertain  disturbances  and  commands  are  modelled  in 
a  realistic  way.  Furthermore,  the  control-input  constraint  and  the  robust  stability  constraints 
for  both  of  the  unstructured  and  parametric  plant  uncertainties  arc  also  taken  into  account 

The  fourth  further  research  suggestion  is  fast  computational  algorithm  for  the  real 
SSV.  Parametric  perturbation  is  a  kind  of  plant  uncertainties  which  occurs  very  often.  The 
closed-loop  system  is  required  to  be  stable  for  all  possible  parametric  perturbations.  It  is 
essential  to  have  a  robust-stability  checking  tool  for  testing  if  the  closed-loop  system  is 
stable  for  all  possible  perturbations  in  a  given  parameter  domain.  The  real  SSV  is  not  only 
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a  robust- stability  checking  tool  but  also  a  stability-margin  measure  which  gives  the 
maximal  allowable  parameter  perturbation  domain  in  which  the  system  remains  stable.  In 
the  design  of  robust  optimal  controller,  we  need  to  compute  the  real  SSV  repeatedly  until 
the  controller  updating  loop  is  terminated.  A  fast  computational  algorithm  for  the  real  SSV 
is  indispensable  in  the  analysis  and  design  of  robust  control  systems. 

For  easy  controller  implementation,  we  would  like  the  order  of  the  controller  to  be 
as  small  as  possible.  The  fifth  further  research  suggestion  is  the  development  of  a  controller 
order  reduction  technique  which  guarantees  the  closed-loop  performance.  Most  of  the 
earing  controller  order  reduction  techniques  are  indirect  methods  in  which  the  basic 
concept  is  to  find  a  reduced-order  controller  which  approximates  the  robust  optimal 
controller  and  then  check  the  closed-loop  characteristics  of  the  reduced-order  controller. 
The  drawback  of  these  indirect  methods  is  that  the  closed-loop  properties  are  not 
considered  during  the  controller  reduction.  The  proposed  approach  is  to  characterize  the  set 
of  all  reduced-order  stabilizing  controllers  and  then  to  find  a  controller  in  this  set  such  that 
the  H°°  cost  function  is  minimized.  This  approach  is  a  direct  method  in  which  the  closed- 
loop  properties  are  always  guaranteed. 
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APPENDIX  A 


Proof  of  Theorems 

Proof  of  Theorem  3.3-2: 

Let  Ng(s)Dg(s)‘1  with  deg  IE>G(s)l  =  n  and  DK(s)‘1NK(s)  with  deg  IDK(s)l  =  n+m 
be  a  right  MFD  (matrix  fraction  description)  of  G22(s)  and  a  left  MFD  of  K(s) 

respectively  [64].  It  is  well  known  that  the  characteristic  polynomial  of  the  closed-loop 
system  is 

^closed-loop^  =  I  Dk(s)Dg(s)  -  NK(s)NG(s)  I  (Al-1) 

Let  L^s)  be  a  greatest  common  left  divisor  of  DK(s)  and  NK(s).  That  is, 

Dk(s)  =  Lk(s)6k(s)  ,  Nk(s)  =  Lk(s)^k(s)  (A  1-2) 

where  £>K(s)  and  &K(s)  are  left  coprime.  It  is  easy  to  see  that  the  zeros  of  Lk(s)  are  the 
uncontrollable  or  the  unobservable  poles  of  the  controller  realization  in  Fig.  3.2-3.  That 
is, 

{  zeros  of  Lk(s)  }  =  ?removal  (Al-3) 

Plug  (A  1-2)  into  (Al-1),  we  have 

fclo«dWs)  =  1  LK<S>  1  • 1  6K(s)DC<s>  -  <VS>N0«  1  CAM) 

The  zeros  of  I  £>K(s)DG(s)  -  ftK(s)NQ(s)  I  are  the  poles  of  the  closed-loop  system  with  a 
minimal  controller  realization.  From  Theorem  3.3-1,  (Al-3)  and  (A  1-4),  we  can  see  that 
(P  +  JP  =  (P  +ip  CA1-51 

closed-loop  with  min.  controller  removal  regulator  observer  Q(s)  v  ' 

To  complete  the  proof,  we  need  to  show  that  ^removal  is  a  subset  of  ^regulator  + 
^observer  w^en  Q(s) ls  minimal.  The  dynamic  equations  of  the  observer- based  controller 
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in  Fig.  3.2-3,  i.e.,  the  block  diagram  inside  the  dotted-line  box,  can  be  written  as 
follows: 


x  =  (A+B2F+HC2+HD22F)  x  + 


[-H  -(B2+HD,2)]| 


F 

-(C2+D22F) 


0  -II  y 


1  D22  U2 


(Al-6a) 


(Al-6b) 


Assume  that  the  added  dynamics  Q(s)  is  described  by  the  following  minimal  realization: 

k  =  A  k  +  B  y  (Al-7a) 


u,  =  Ck  +  Dy 


(Al-7b) 


The  controller  K(s)  is  just  a  combination  of  (A  1-6)  and  (A  1-7).  From  (A  1-6)  and  (Al- 
7),  we  have  the  dynamic  equations  of  the  controller  K(s)  as  follows: 


where 


all 

«12 

a21 

tt22 

Ta 

u  =  [y,  Y2]  +  8y 


Pj  =  -  H  -  (B2+HD22)  a-DD^-'D 


P2  =  B  +  B  D22(I-DD22)  D 


Yj  =  F  +  (Cj+D^F) 


(Al-8a) 


(Al-8b) 


(Al-9a) 

(Al-9b) 

(Al-9c) 
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Y2  =  -  (I-DD^C 

(Al-9d) 

an  =  A  +  HC2+(B2+HD22)y, 

(Al-9e) 

=  A  +  B2F  -  pj  (C2+D22F) 

(Al-9f) 

ai2  =  (B2+HD22)  ^2 

(Al-9g) 

°2i  =  -P2(C2+D22F> 

(Al-9h) 

a22  =  ^  ®  ^22  ^2 

(Al-9i) 

8  =  -a-DD^D 

(A1-9J) 

Now,  assume  that  the  state-space  representation  (A  1-8)  of  the  controller  K(s)  is 
unobservable.  Then  by  PBH  test  [64],  there  exists  a  nonzero  vector  t,  such  that 


Si' 

^2 


an  ai2 

V 

V 

=  X 

_  tt21  tt22  _ 

^2 

^2_ 

=  * 


(Al-10a) 


(Al-lOb) 


for  some  eigenvalue  X  of  (Al-8).  Note  that  it  is  the  eigenvalue  X  that  is  unobservable. 
From  (Al-lOb),  we  get 


“n^i  +  ai2^2  =  ^  .  (Al-lla) 


which,  by  using  (Al-9e)  and  (Al-9g),  is  rewritten  as 


(A  +  HC2)^  +  (B2+HD22)(y1^1+Y2^2)  =  (Al-llb) 

In  view  of  (Al-lOa),  the  above  equation  reduces  to 
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( A  +  HC2)^  =  x^ 


(Al-llc) 


which  clearly  establishes  that  the  unobservable  eigenvalue  belongs  to  ^obMiVW- 


Proceeding  similarly,  it  can  be  shown  that  if  (A  1-8)  is  uncontrollable  then  the 
uncontrollable  eigenvalue  belongs  to  ^regulator  • 


Note  in  the  above  development  that  ^  =  0  contradicts  the  minimality  assumption  of 
Q(s).  Thus, 


V 

V 


unobservable 

uncontrollable 


C 

C 


V 


observer 

V 

regulator 


(Al-12a) 

(Al-12b) 


where  IP 


unobservable 


is  the  set  of  all  the  unobservable  poles  of  the  controller  K(s). 


<P 


uncontrollable 


is  defined  similarly.  This  completes  the  proof  of  Theorem  3.3-2. 


Proof  of  Remark  4.2-1: 

Clearly,  the  realization  (4.2-4)  is  controllable  if  the  pair 

{(A2  -  L2L2X2)T,  (D±W2  -  D12Rd1B2U2X2)T)  (A2-la) 

is  controllable.  This  pair  is  controllable  if  there  exists  a  positive  definite  &  such  that 

<A2  -  L^x/x  +  X  (Aj  -  LjL^Xj) 

+  (D  W2-D,2Rd'b^U2X2)T  (D  Wj-D.jRo'b^UjXj).  o  (A2-lb) 

Substituting  by  X2,  then  the  left  hand  side  of  (A2- lb)  can  be  reduced  to 

A2X2  +  X2A2  -  X2L2L2X2  +  W2W2  (A2-lc) 

which  is  zero  according  to  (4.2-3e). 
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Proof  of  Theorem  4.3-7: 


Recall  that 


M-(Y)  ••= 


inf  ll(R  +Q)(fl-R*  R  )'1/2I 

£>e  (RH°°)mxr 


(A3-1) 


The  theorem  will  be  proved  by  a  contradiction.  Let 


a  =  V  H2(Y)  Y2  +  [  1  -  U2(Y)  1  c2 


(A3-2) 


Suppose  that  (4.3-14)  is  not  true.  Then  from  Theorem  4.3-1,  Y0  <  a  implies  that 

H(a)  :=  inf  II  (Rn+$)  (a2I-R*R21  )'1/2II  _  <  1  (A3-3) 

Qe  (RH°°)mxr  2  21 

An  optimal  Q  for  (A3-3)  is  denoted  by  i.e., 

P(a)  =  ll(Rn+01)(a2I-R*21R2iy1/2lloo  <  1  (A3-4) 


Thus, 


(a2I-R21R21  )1/2*  (R„+0,)*  (Rn+^)  (a2I-R21R21  ym  <  I  for  all  ©  (A3-4) 

<=>  (Rn+ftj)*  (R11+01)  <  a2I-R;,R21 

=  lt2(Y)  Y2!  +  D-H2(Y )1  c2 1  -  R2iR2i 
£  H2(Y)  Y2  I  +  [1-M-2(Y>]  R2lR21  -  R21R21 
=  M-2(Y)  ( Y2 1  -  R2iR2i )  for  all  o)  (A3-5) 

»  (Y2!-R21R21  )~,/2*(ri  ,+01)*(r1  1  +0i)<T^I-R21R2i  )'1/Z  <  |A2(Y)I  for  all©  (A3-6) 
**  ll(Rn+^i)(Y2I-R2iR2i)'1/2HM  <  m  (A3-7) 


which  contradicts  the  fact  that  |i(y )  is  the  infimum  as  stated  by  (A3-2).  Hence, 


Y°  *  a  ~  V  4%  Y2  +  [  1  -  M-2(Y)  1  c 


.2,  .  „  2 


(A3-8) 
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Proof  of  Corollary  4.3  8: 


Denote  the  right-hand  side  of  (4.3-15)  (  or  (4.3-17) )  by  p3.  The  function  |i(y)  can  be 
considered  as  a  curve  in  a  two-dimensional  space.  The  point  ( y  ,  1  )  is  the  intersection 

of  the  curve  y  =  |i(Y)  and  the  horizontal  line  y  =  1.  Draw  a  straight  line  connecting  the 
two  points  (  Yr  M.(Yj) )  and  (  y5,  |i(Y5)  )•  This  straight  line  will  cross  the  horizontal 
line  y  =  1  at  the  point  (P3,  1).  Since  the  function  p(y)  is  continuous,  convex,  and 

strictly  monotonically  decreasing,  it  is  easy  to  see  that 

Y0  >  P3  if  H(Y5)  >  1  and  yo  <  P3  if  ^(Y5)  <  1- 
Inequalities  (4.3-16)  and  (4.3-18)  can  be  proved  in  a  similar  way. 


Proof  of  Theorem  5.2-3: 

Without  loss  of  generality,  we  will  show  that  a.'s,  i=0,l,2,...,n,  are  linear  (affine) 
functions  of  8,  while  82,  ...,  8m  remain  unchanged.  Let  m;.(s)  be  the  i-j  entry  of  the 

matrix  M(s).  Then 

det  [  I  +  M(s)A  ] 

1+mjjSj  mi2^2  •  •  *  mim^m 

102,8,  l+m2282  •  •  • 

mml^l  mm2^2  '  '  '  i+mmnAn 


1 

m,2S2 

■  ■  •  “im8" 

m,. 

m1282  .  .  . 

m.m8m 

0 

1+m22®2 

■  •  ■  “to8" 

+  Si 

1+m2282  .  . 

0 

mm282 

■  ■  l+«W8m 

mmA  •  •  • 

:=  foj(s)  +  8j  ^(s) 

Here  frij (s)  that  fh^s)  are  rational  functions  with  real  coefficients.  Now,  we  can  write 

ftv(s)  +  6,  «Us)  =  n'(S)  +  8,  ■  ,  "  Ws) : 

1  2  dc(s)dx(s)  01  d^s^s)  dc(s)d1(s)d2(s) 

where  the  (n^s)  ^(sidjfs))  and  {n2(s),  dc(s)d2(s) }  are  coprime  polynomial  pairs,  and 
dc(s)  is  a  greatest  common  factor  of  the  denominator  polynomials  of  fhj(s)  and  fo^s). 
The  polynomials  n^s^fs)  and  n2(s)d1(s)  can  be  represented  as 

aQsn  f  ajS"'1  +  a2sn'2  + . +  an 

and 

b„sn  rb,sn  l  +  b,sn'2  + . +  b 

0  12  n 

respectively.  Therefore,  the  characteristic  polynomial  is 
aQsn  +  ctjs"'1  +  a2sn"2  + . +  an 

=  sn  +  (aj+Sjbj)  sn  l  +  (a^Sjbj)  s"'2  + . +  (an+8jbn) 

i.e., 

a.  =  a;  + 8jb.  i  =  0,l,2,...,n 

Proof  of  Theorem  5.2>6: 

Recall  that  3  =  {  8  :  -1  <  8;  ^  1,  i  =  1,2,... ,m  }.  To  complete  the  proof  we  need  to 
show  that  for  any  8*  e  3,  the  image  of  5*  is  inside  the  polytope  IP (3)  if  the  mapping 
is  multilinear.  Let 

8  =  [8,  ,  82  ,  83  ,  ....  8m  ]T  e  3 

and  define 

5  00  =  [8,  ,  82  ,  ....  ^  ,  ak+1,  ak+2 . am]T 
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where  0  *"  k  <,  m,  and  <x  is  either  +!  or  -1,  for  k+1  <,  i  <,  m.  Note  that  8*(k)  represents 
one  of  the  2m'lc  points  which  associate  with  8*  on  the  boundary  hyperplanes  of  J§.  8*(0) 
is  one  of  the  vertices  of  J§  and  8*(m)  =  8*.  We  shall  prove  by  mathematical  induction 
that  ;J(8*(m))  e  1P(J§). 

Assume  that  the  image  of  the  2m  k  S*(k)'s  are  all  in  1P(JS))  and  define  * 

6  (k  )  =  [8,  ,  82  ,  \  .  +1,  ak+2’  °rJ 

and 

8*(k  )  =  [8j*,  82*,  Sk*.  -1,  ak+2, ....  crm]T 

Let  the  line  segment  connecting  8*(k+)  and  8*(k  )  be  denoted  by  ik.  Obviously,  ik  is 
parallel  to  8k+1-axis,  since  the  only  coordinate  change  from  8*(k+)  to  8*(k  )  is  along  the 
8k+1-axis.  Therefore,  the  image  of  is  also  in  P(J§)  by  the  virtue  of  Lemma  2.5.  Since 
8* (k+1)  is  a  point  on  ik,  its  image  must  also  in  1P(JB).  The  fact  that  the  image  of  8*(0) 
is  in  iP(jB)  is  established  by  the  definition  of  <P(  £)).  This  completes  the  proof. 

Proof  of  Theorem  5.2-7: 

The  proof  is  similar  to  De  Gaston  and  Safonov’s  [48].  The  only  difference  is  that  their 
co-domain  is  the  complex  plane  while  ours  is  tl.  c  coefficient  space.  The  roots  of  the 
characteristic  equation,  i.e.,  the  eigenvalues  of  A-BAC  are  uniquely  determined  by  the 
parameter  vector  8  and  therefore  so  is  the  coefficient  vector  a.  For  any  8*  e  JS)V  there 
exists  a  8y  e  jj)2  such,  that  8y  =  8X  and  thus  0(8y)  =  $(8*).  Hence  c  SJ(JS2). 

Since  ^(JSj)  is  the  smallest  convex  set  that  contains  all  the  points  of  d(J®  j),  (5.2-6)  is  a  ‘ 

direct  consequence  of  (5.2-5). 
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Proof  of  Theorem  5.2-8: 


Again,  the  proof  is  similar  to  De  Gaston  and  Safonov's  [48].  The  vertices  of  £)x  and  JS)2 
are  mapped  into  P(J§)  and  therefore  P(JS)j),  V(js>2),  and  P(J§j)  vj  1P(^2)  are 
contained  in  P(JD).  Since  d(j0j)  <=  JP^),  $(j§2)  c  F(J§2),  and  u  d(J§2)  = 

d(J9),  we  have  d(j®)  c  u  1P(«©2). 
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APPENDIX  B 


Real  Solution  Existence  Checking  by  Using  Sturm  Theorem 
Given  a  polynomial  Pn(x)  of  degree  n  and  an  interval  [a,b]  on  the  real  axis,  find 


how  many  roots  of  Pn(x)  are 

located  inside  [a,b]  ? 

Sturm  Sequence: 

Let 

and 

Pn(x)  =  anxn  +  a.x"'1  +  ...+  a  .  x  +  a 
n'  '  0  1  n-1  n 

(B-l) 

V — 

II 

'ST 

'w' 

(B-2) 

(B-3) 

Divide  fQ  by  fj  and  denote  the  remainder  by  -fr  Divide  fj  by  f2  and  denote  the  remainder 
by  -f3,  continue  the  process  until  fn  is  a  constant  To  illustrate, 

f0  =  Ql  f 2 
fl=Q2f2-f3 

(B-4) 

*in-2  =  9n-l  *n-l  *  *n 

fn-l=Q„f„ 

To  actually  compute  the  coefficients  of  the  polynomials  f.,  we  suggest  forming  the 

following  table:  ♦ 
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where  the  a.'s  (i=0,l,2,...,n)  are  the  coefficients  of  the  original  polynomial  fQ(x).  The  b.'s 
(i=0,l,2,...,n-l)  are  the  coefficients  of  ft(x).  The  c/s  (i=0,l,2,...,n-l)  are  given  as 

=  boai+i "  ^i+i  (B-5) 

and  the  c.'s  are  the  coefficients  of  f2(x)  which  are  given  as 

ci =  5obi+i '  b<A+i  (B-6) 

Similarly  the  coefficients  of  f.(x)  (i  =  3,  4, ....  n)  can  be  computed.  Once  these  functions 
are  constructed,  it  can  be  shown  that  they  form  a  Sturm  sequence  [63]  with  the  following 
properties: 

i)  For  any  value  xQ  e  [a,b]  at  which  f.(x)  =  0,  we  have  fi.1(x0)-fi+l(x())  <  0. 

ii)  The  last  function  in  the  sequence  fn(x)  does  not  vanish  for  any  x  e  [a,b]  .  Also, 
it  is  required  that  fQ(x)  does  not  vanish  at  the  end  points  of  the  interval  [a,b]. 

Note  that  any  given  row  of  the  above  table  can  be  multiplied  or  divided  by  any  nonzero  real 
constant  without  affecting  the  properties  of  the  functions  f.(x). 

Theorem  (Sturm): 

Let  the  functions  f0(x),  fj(x),  f2(x),  ...,  fn(x)  be  a  Sturm  sequence  and  let  V(x) 

represent  the  number  of  variation  of  sign  of  the  sequence  at  x,  then  the  number  of  real  roots 
in  the  interval  [a,b]  is  exactly  V(a)  -  V(b)  with  repeated  roots  counted  only  once. 

A  proof  of  this  theorem  can  be  found  in  [62]. 
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